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ABSTRACT: 

Steady  shock-free  transonic  compressible  fluid  flow  about 
a  series  of  symmetric  airfoils  without  circulation  is  found  by 
computing  analytic  solutions  of  the  hodograph  equations  for  an 
ideal  polytropic  gas  in  two  dimensions. 

The  hodograph  variables  are  extended  analytically  into 
complex  four-space  and  the  equations  are  written  in  character- 
istic normal  form.   This  procedure  has  been  used  previously  on 
the  detached  shock  problem  by  P.  Garabedian  [7],  E.  Swenson  [24], 
and  D.  Korn  [12].   Characteristic  coordinates,  s  and  t,  analogous 
to  w  =  u  -  iv  and  w*  =  u  +  iv  for  Laplace's  equation 


x   „  =  0 

WW* 


governing  incompressible  flow  are  found  and  they  are  used  to 
generalize  the  singular  solution  of  the  form 

x  =  Re 


/T=w 

for  incompressible  flow  around  an  airfoil  with  unit  speed  at 
infinity.   The  analogous  singular  solution  for  compressible  flow 
has  the  form 


Re 


A(t) 


+  B(s,t) 


where  A  and  B  are  regular.   A  large  class  of  airfoils  can  be  found 
by  adding  a  regular  solution  of  the  characteristic  initial  value 
problem  to  this  singular  solution.   The  problem  of  continuing 
this  solution  into  the  supersonic  region  is  achieved  by  choosing 
paths  on  the  initial  characteristics  such  that  the  complex  solu- 
tion winds  around  the  complex  sonic  line  (the  set  of  points  in 
the  complex  hodograph  where  the  speed  q  is  equal  to  the  critical 
speed  c#)  and  end  up  at  real  points  in  the  hodograph.   In  fact, 
one  path  can  be  found  which  will  yield  the  solution  in  the  whole 
supersonic  region  outside  the  body. 

The  method  of  Massau,  used  previously  to  obtain  solutions 
of  the  detached  shock  problem,  is  used  to  solve  a  characteristic 
initial  value  problem.   This  second  order  scheme  integrates  the 
equations  along  complex  characteristics  and  hence  stability  is 
assured.   The  equations  become  singular  at  the  sonic  line  and 
hence  there  is  a  loss  of  accuracy  there.   Accuracy  is  checked  by 
taking  a  deferred  approach  to  the  limit.   In  general,  a  solution 
can  be  computed  in  a  matter  of  four  or  five  minutes  on  the 
CDC  6600.   Points  on  the  body,  the  pressure  distribution  along 
the  body,  as  well  as  the  characteristics  in  the  supersonic  zone, 
are  plotted  on  a  Calcomp  plotter.   Solutions  which  result  by 
varying  a  few  parameters  in  the  characteristic  initial  data  are 
found  quickly  with  the  aid  of  a  time-sharing  system. 

By  varying  the  data  given  along  the  characteristics,  we 
are  able  to  control  the  shape  of  the  nose  and  tail  rather  easily. 
In  particular,  inserting  a  logarithmic  singularity  near  the  nose 
or  tail  is  an  effective  way  to  change  the  flow  chiefly  in  the 


VI 


neighborhood  of  the  singularity.   Similarly,  the  location  and 
occurrence  of  limiting  lines  can  be  controlled.   By  placing  a 
limiting  line  just  inside  the  body  near  the  front  we  can  get  a 
"peaky"  pressure  distribution  as  described  by  H.  Pearcey  [20], 
Other  types  of  pressure  distributions,  perhaps  with  inferior 
aerodynamic  qualities  are  also  obtained. 


Vll 


1.  Introduction 

In  this  paper  we  investigate  a  technique  of  computing 
airfoil  sections  that  have  shock-free  transonic  flow  around 
them.   By  transonic  we  mean  that  the  speed  of  the  aircraft  is 
less  than  the  speed  of  sound,  but  close  enough  to  it  so  that 
on  top  of  the  wing,  where  the  airflow  is  fastest,  the  speed  is 
greater  than  the  speed  of  sound.   For  air  the  Reynold's  number 
is  extremely  high,  and  hence  the  viscous  effects  will  be  con- 
fined chiefly  to  a  thin  boundary  layer  provided  separation 
effects  can  be  suppressed.   We  neglect  the  viscous  terms  and 
take  for  our  equations  the  equations  of  potential  flow. 

Thus,  the  problem  is  to  find  smooth  transonic  steady- 
state  solutions  to  the  equations  of  irrotational  motion  of  an 
inviscid,  compressible  fluid  about  an  airfoil  section  in  two 
dimensions.   Two  main  difficulties  arise  in  analyzing  this 
problem.   First  of  all,  the  partial  differential  equations  of 
motion  are  nonlinear.   This  difficulty  can  be  avoided  by  using 
the  hodograph  transformation,  since  the  equations  become 
linear  in  the  hodograph  plane.   Secondly,  when  the  speed  be- 
comes supersonic,  the  equations  change  type,  from  elliptic  to 
hyperbolic.   It  has  been  shown  by  Morawetz  [l6]  that  smooth 
transonic  solutions  to  the  equations  of  motion  do  not  exist  for 
all  airfoil  shapes.   It  follows  that  the  problem  of  computing 
shock-free  transonic  flow  past  a  given  profile  is  not  well  posed. 
We  overcome  this  second  difficulty  by  solving  an  inverse  problem. 
That  is,  we  compute  a  smooth  transonic  flow  and  then  find  the 


body  which  generates  this  solution.   To  find  a  smooth  solution 
we  extend  all  the  variables  into  the  complex  plane,  as 
suggested  by  Garabedian  [8]  and  used  successfully  on  the  de- 
tached shock  problem  (cf.  [7],  [12],  (24]).   Then  we  solve  a 
characteristic  initial  value  problem  along  the  complex  charac- 
teristics.  This  procedure  can  be  carried  out  on  an  electronic 
computer  by  a  simple  numerical  scheme. 

In  view  of  recent  experimental  success  (cf.  [11],  [20], 
[25])  in  achieving  virtually  shock-free  flow  in  wind  tunnel 
tests,  the  problem  has  gained  new  impetus.   Function-theoretic 
methods  based  on  the  hodograph  transform  were  first  used  by 
Chaplygin  in  1902  in  his  work  on  gas  jets.   Such  methods  have 
been  investigated  by  other  authors.   Von  Karman  and  Tsien  [10] 
and  Tomotika  and  Tamada  [25]  have  constructed  solutions  by 
modifying  the  equation  of  state.   Also,  solutions  have  been 
found  by  suitably  superposing  particular  solutions  of  the  hodo- 
graph equations  which  are  found  by  separation  of  variables. 
Series  solutions  with  coefficients  related  to  an  incompressible 
flow  are  used  and  hence  the  problem  of  convergence  and  analytic 
continuation  of  these  series  is  involved.   Goldstein,  Lighthill, 
and  Craggs  [9]*  and  Cherry  [3]  were  the  first  to  develop 
successful  methods  for  the  analytic  continuation  of  the  solu- 
tion for  the  case  related  to  non-circulatory  incompressible  flow 
around  a  circular  cylinder.   Later,  Lighthill  [14]  presented  a 
more  general  form  of  the  theory  in  which  he  defines  an  operator 
which  assigns  a  compressible  flow  to  each  incompressible  one. 


Bergman  [l]  has  presented  a  theory  based  on  his  general 
method  of  mapping  analytic  functions  into  the  space  of  solu- 
tions of  the  hodograph  equations.   He  uses  the  analytic  proper- 
ties of  the  equations  to  construct  the  fundamental  solution, 
which  is  used  as  the  kernel  of  his  integral  operator.   Numeri- 
cal solutions  can  be  found  by  expanding  the  integrand  as  a 
series  and  performing  the  integration.   In  some  ways  his  method 
resembles  ours.   In  fact,  his  idea  of  determining  the  funda- 
mental solution  has  been  incorporated  in  the  method  used  in 
this  paper. 

Perhaps,  the  most  recent  work  on  this  problem  employing 
the  hodograph  method  is  the  work  done  at  the  National  Aerospace 
Laboratory  (NLR)  in  Amsterdam.   Nieuwland  [17]  has  expanded  on 
the  work  of  Lighthill  and  has  developed  a  technique  of  computing 
a  three-parameter  family  of  thin  quasi-elliptical  airfoil 
sections.   He  pieces  together  series  of  solutions  for  the  stream 
function  found  by  separation  of  variables.   He  then  computes 
the  velocity  potential  and  the  physical  coordinates  x  and  y. 
With  the  aid  of  a  one  parameter  family  of  solutions,  he  is  able 
to  ensure  that  x  and  y  are  single  valued.   A  catalogue  of  non- 
circulatory  flows  arrived  at  by  this  method  is  presented  by 
Boerstoel  [2].   Experiments  in  a  wind  tunnel  with  these  solu- 
tions by  Spee  and  Uijlenhoet  [25]  show  agreement  with  the  theory, 
thus  proving  the  validity  of  these  solutions.   They  show  further 
that  the  solutions  obtained  in  this  manner  have  a  certain  amount 
of  stability  to  fluctuations  in  speed  and  angle  of  attack. 


The  method  used  in  this  paper  is  similar  to  the  above 
techniques  in  that  they  all  use  an  inverse  method.   Similarly, 
the  equations  are  made  linear  by  using  a  hodograph  transforma- 
tion.  In  other  ways  the  methods  differ  widely.   Our  method 
uses  finite  differences  to  solve  a  characteristic  initial  value 
problem  in  complex  space  as  opposed  to  computing  series.   We 
will  see  that  the  whole  supersonic  region  can  be  computed  by 
choosing  one  set  of  initial  paths.   The  characteristics  of  the 
flow  are  readily  obtained  by  our  procedure  and  thus  the  occur- 
rence and  location  of  limiting  lines  is  easily  detected.   We 
use  x  and  y  as  the  dependent  variables  because  we  do  not  rely 
on  the  use  of  separation  of  variables.   By  doing  this,  we  avoid 
computing  the  velocity  potential.   The  single-valuedness  of  the 
solution  in  the  physical  plane  is  assured  by  the  construction 
of  the  singular  solution.   With  our  method  we  are  able  to  com- 
pute solutions  for  any  set  of  initial  data  in  about  5  minutes  of 
machine  time.   The  solutions  obtained  have  a  more  favorable 
thickness  to  chord  ratio  for  a  given  Mach  number  than  the  ones 
obtained  by  the  other  methods. 

In  Section  2  we  write  down  the  equations  of  motion  and 
we  introduce  the  hodograph  transformation.   All  of  the  variables 
are  extended  into  the  complex  domain  in  Section  ~$.      Character- 
istic coordinates  are  introduced  in  Section  4  and  the  equations 
of  motion  are  reduced  to  characteristic  normal  form.   Section  5 
is  a  digression  on  incompressible  flow  about  an  airfoil.   The 
problem  of  incompressible  flow  is  formulated  and  solved  in 
terms  of  a  characteristic  initial  value  problem  in  the  hodograph 


plane  which  has  been  extended  by  analytic  continuation.   The 
results  of  this  section  motivate  much  of  the  work  for  the 
compressible  case.   In  Section  6  we  show  how  to  find  an  airflow 
by  combining  a  regular  solution  and  a  singular  solution  which 
is  related  to  the  fundamental  solution.   A  particular  choice  of 
characteristic  coordinates  is  made  in  analogy  with  the  formula- 
tion in  Section  5.   A  closed  form  solution  of  the  characteristic 
transformation  is  presented.   Simplifications  which  result  from 
neglecting  circulation  are  discussed.   Section  7  contains  a 
description  of  the  finite  difference  scheme  used  to  compute  solu- 
tions of  the  equations  on  an  electronic  computer.   We  show  how 
to  choose  paths  of  initial  data  which  yield  the  solution  in  the 
subsonic  portion  of  the  flow.   We  introduce  the  so-called  sonic 
locus  (cf.  [25])  in  Section  8  and  explain  how  it  is  used  to 
continue  the  solution  into  the  supersonic  zone.   We  show  how  to 
choose  one  set  of  initial  characteristic  paths  which  yield  the 
whole  supersonic  zone  outside  the  body.   A  brief  discussion  of 
the  singularities  that  can  arise  in  the  hodograph  transformation 
appears  in  Section  9,    and  we  explain  how  these  singularities  may 
be  used  to  form  a  cusp  at  the  trailing  edge  of  the  airfoil  and 
a  "peaky"  pressure  distribution  as  described  by  Pearcey  [20]. 
The  numerical  results  obtained  on  a  CDC  6600  for  the  case  without 
circulation  are  presented  in  Section  10.   Accuracy  as  well  as 
computation  time  and  storage  requirements  are  discussed  here. 
Many  examples  of  flows  computed  by  this  method  are  illustrated. 
Section  11  discusses  the  possibility  of  physically  realizing 


flows  obtained  by  this  method.   Finally,  a  copy  of  the  coded 
Fortran  program  is  included  in  the  back  of  this  paper. 

I  would  like  to  express  my  appreciation  to  Professor 
Paul  R.  Garabedian  who  suggested  this  problem.   His  enthusiasm, 
advice,  and  patience  proved  invaluable  to  me .   I  would  also 
like  to  thank  H.  Samoraj  for  typing  this  paper  and  A.  Weiner 
who  prepared  an  excellent  set  of  pen  and  ink  drawings  for  the 
figures. 


2.  The  Equations  of  Motion  and  the  Hodograph  Transformation 

The  equations  of  motion  of  an  inviscid,  isentropic,  poly- 
tropic  gas  for  steady,  irrotational  flow  in  two  dimensions  are 
presented  by  several  authors  [4,  8,  15] •   To  be  consistent  with 
the  literature  we  let  x  and  y  be  Cartesian  coordinates  and  let 
u  and  v  be  the  velocity  components  in  the  corresponding  direc- 
tions.  The  density  will  be  denoted  by  p  and  the  pressure  by  p. 
The  equation  of  state  for  a  polytropic  gas  is 


(2.1)  p  =  Ap 
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where  A  and  y   are  constants.   For  air  y  =   1.4.   The  local  sound 
speed  is  defined  as  (-p)  '      and  is  denoted  by  c.   With  this 
notation  the  equations  of  motion  are 


(2.2a)        (c2-  u2)u  -  2uvu  +  (c2  -  v2)v   =  0 


(2.2b)  vx  -Uy  =  0 


We   also  have  Bernoulli's   equation 


,_,,  222        „2  2  2,2 

(2.3)  q     +yzi    c=q;         q=u+v 


o 
where  q  is  the  limit  speed  so  that  c  ,  p,  and  hence  p  are  func- 

2 
tions  only  of  q  . 

The  fact  that  mass  is  conserved  is  expressed  by  the  con- 
tinuity equation 
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(2.4)  (PU)   +  (PVKr    =     ° 

-A-  J 

which  has  been  used  in  the  derivation  of  (2.2a).   This  equation 
implies  the  existence  of  a  stream  function  f   such  that 


(2-5)  PQ^X  -  pv  ;    pQfy   =  -pu 


where  p   is  the  stagnation  density. 

The  coefficients  of  the  homogeneous  quasi-linear  system 
(2.2a-b)  are  functions  of  only  u  and  v  and  hence  it  is  possible 
to  transform  this  system  into  an  equivalent  linear  system  by 
interchanging  the  roles  of  the  dependent  and  independent  vari- 
ables.  In  any  region  where  the  Jacobian 


(2.6)  j=uv-uv 

x  y   y  x 


is  non-zero  we  may  consider  x  and  y  as  functions  of  u  and  v  for 
a  solution  of  (2.2a-b).   Using  standard  formulas  of  advanced 
calculus  we  obtain  the  relations 


(2.7a)  ux  =  jyy  ;   uy  =  -jxy 


(2.7b)  v  =  -jy   :    v  =   jx 

'  x    °    v      y      u 


Inserting  these  relations  into  (2.2a-b)  we  obtain  the  linear 
system 

(2.8a)  (c2  -u2)yy  + 2uvxy  +  (c2-  v2)xu  =  0 

(2.8b)  _yu  +Xy   =   0 
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Similarly  if  the  Jacobian 


(2'9)  J  =  Xuyv-Xvyu 


does  not  vanish,  then  any  solution  of  (2.8a-b)  is  also  a  solu- 
tion of  (2.2a-b) . 

Any  transformation  which  interchanges  the  dependent  and 
independent  variables  is  called  a  hodograph  transformation. 
Note  that  solutions  for  which  j  =  0  cannot  be  found  in  this  way. 
Included  in  this  category  are  the  simple  waves.   Solutions  with 
limiting  lines  on  which  J  =  0  will  be  discussed  in  Section  9. 


3-  Complex  Extension 

System  (2.8a-b)  can  be  written  in  the  matrix  notation 


which  is  of  the  form 


(3.2) 


SX  +TX   =  0 
v    u 


The  type  of  such  a  second  order  system  is  determined  by  examining 
the  two  roots,  A  and  A  ,  of  the  equation 


(3-3) 


Det(SA  -  T)  =  0 


from  which  we  find  that 


(3-4) 


*+  = 


,  pk      2 

uv  ±  c/q  -  c 


2    2 

c  -  u 


When  both  roots  are  real  and  unequal  the  system  is  hyperbolic, 
while  if  the  two  roots  are  complex  conjugates  the  system  is 
elliptic. 

We  call  the  dimensionless  quantity  —  the  Mach  number  and 
denote  it  by  M.   Clearly  both  roots  of  (3-3)  are  real  and  unequal 
when  M  >  1,  while  for  M  <   1  the  two  roots  are  complex  conjugates. 
For  M  =  1  the  two  roots  are  real  and  equal.   From  Bernoulli's 
equation  (2.3)  we  see  that  the  locus  of  points  in  the  hodograph 
plane  where  M  =  1  is  a  circle,  q  =  constant.   We  call  this 


10 


constant  the  critical  speed  and  denote  it  by  c*.   The  critical 
speed  is  related  to  the  limit  speed  by  equation 


(3-5)  c„  =  nq  ;    u.' 


7-1 
7+1 


For  q  >   c^,  M  >•  1  and  the  flow  is  said  to  be  supersonic,  while 
for  q  <   c^,   M  <   1,  the  speed  is  said  to  be  subsonic. 

For  supersonic  flow  A   and  A_  are  real  and  the  equations 
are  hyperbolic.   In  this  case  the  initial  value  problem  (x  and 
y  assigned  along  any  non-characteristic  curve)  and  the  charac- 
teristic Initial  value  problem  (x  or  y  assigned  on  one  charac- 
teristic of  each  family)  are  well  posed.   For  subsonic  flow  A 
and  A_  are  complex  conjugates  and  the  equations  are  elliptic. 
The  initial  value  problem  and  the  characteristic  initial  value 
problem  are  no  longer  well  posed  In  this  case,  and  we  generally 
assign  boundary  data  to  formulate  a  properly  posed  problem. 

We  are  concerned  here  with  transonic  flow,  that  is,  flow 
which  is  partly  subsonic  and  partly  supersonic.   Therefore,  it 
is  not  immediately  obvious  what  type  of  initial,  characteristic 
initial,  or  boundary  conditions  will  comprise  a  well  posed 
problem.   A  correct  problem  has  been  formulated  by  Tricomi  [27] 
and  Frank '1  [6]  by  studying  the  nature  of  a  somewhat  simpler 
mixed  equation.   Their  results  include  two  possibilities  of  a 
well  posed  problem  for  an  equation  of  mixed  type. 

There  is  another  alternative.  We  note  that  the  coeffi- 
cients of  system  (3.1)  are  all  analytic  functions  of  u  and  v. 
Using  the  method  of  Garabedian  [8]  we  extend  both  u  and  v  Into 
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the  complex  domain.   In  this  way  the  independent  variables 
range  over  a  four  dimensional  space  instead  of  the  plane. 
Initial  data  or  characteristic  initial  data  may  be  uniquely 
extended  into  this  complex  space  providing  only  that  they  are 
analytic  functions  of  u  and  v.   This  method  of  extension  has 
been  used  in  fluid  dynamics  previously  by  Garabedian  [7], 
Swenson  [24],  and  Korn  [12].   The  existence  of  a  solution  within 
the  class  of  analytic  functions  is  ensured  at  least  locally  by 
the  Cauchy-Kowalewski  theorem.   How  far  we  can  extend  such  a 
solution  will  become  self-evident  after  the  numerical  technique 
is  presented  in  Section  7-   In  complex  space  the  notion  of  type 
plays  no  importance.   Instead,  we  note  from  equation  (3.4)  that 
the  roots  are  distinct  provided  that  q  ^  c^.   The  locus  of 
points  where  q  =  c^  consists  of  a  two  dimensional  manifold  in 
the  four  dimensional  (u,v)-space  which  we  shall  call  the  complex 
sonic  line.   We  shall  see  in  Section  8  that  the  complex  sonic 
line  plays  an  important  role  in  the  continuation  of  our  flow 
into  the  supersonic  region. 

By  formulating  our  problem  as  an  initial  value  problem  or 
a  characteristic  initial  value  problem  we  are  solving  an  inverse 
problem.   That  is,  we  do  not  prescribe  the  shape  of  the  airfoil 
and  compute  the  flow,  rather,  we  compute  an  analytic  solution 
and  find  the  body  which  generates  this  flow.   In  fact,  Morawetz 
[l6]  has  shown  that  it  is  not  possible  to  find  a  smooth  solution 
for  an  arbitrary  airfoil.   Our  method  enables  us  to  find  only 
those  airfoils  for  which  the  airflow  is  smooth. 
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The  real  and  imaginary  parts  of  x  and  y  are  separately 
solutions  if  x  and  y  are  complex  solutions  since  system  (3.1) 
is  linear  and  the  coefficients  are  real  in  the  real  hodograph, 
In  what  follows  we  distinguish  between  the  real  variables  and 
the  extended  complex  variables  only  when  their  meaning  is  not 
clear  from  the  context. 
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4.  Characteristic  Coordinates 

In  the  previous  section  we  expressed  our  system  of  equa- 
tions (J.l)  in  matrix  form  (3.2).  Now  let  x  be  a  non-zero  row 
vector  such  that 

(4.1)  x(SA  -T)  =  0 
and  hence 

(4.2)  xSA  =  xT 

Clearly  a  non-trivial  x  exists  only  if  SA  -  T  is  singular,  i.e. 

(4.3)  Det(SA  -  T)  =  0 


Equation  (3.4)  gives  the  two  solutions,  A  and  A  ,  of  this  equa- 
tion. Multiplying  (3.2)  on  the  left  by  the  corresponding  x  and 
x  and  using  (4.2)  we  arrive  at 


(4.4a)  x+S(Xy  +  A+Xu)  =  0 


(4.4b)  X_S(XV  +A-XJ  =  ° 


which  are  two  scalar  equations  each  of  which  involves  differen- 
tiation in  only  one  direction.   Carrying  out  the  calculations 
we  find  that 

(4.5)  x±S  =  (A±,l) 

We  introduce  the  coordinates  £,  and  r|  with  the  properties 
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v  d£    dv     -  ou 

(^•6b)  tJt  =  ^  +  a+It 

v  or)   av    +  OU 


It  is  easy  to  see  that 


(4.7a)  u^  =  A+v; 


(4.7b)  u^  =  A_v^ 

Using  (4.5)  in  equation  (4.4a-b)  we  get  two  additional  equations 

(4.8a)  y^  +  A_x^  =  0 

(4.8b)  y^  +  A+x^  =  0 

The  coordinates  i   and  r\   are  called  characteristic  coordi- 
nates and  equations  (4.7a-b),  (4.8a-b)  are  the  equations  of 
motion  written  in  terms  of  the  characteristic  coordinates.   The 
curves  |  =  constant  and  r\    =  constant  are  called  characteristics 
and  A  and  A_  are  the  characteristic  directions.   Since  we  have 
extended  u  and  v  Into  the  complex  domain,  it  is  clear  that  i   and 
T)  must  also  be  complex.   Thus  the  characteristics  are  two 
dimensional  surfaces  in  four  dimensional  space. 

Clearly  any  function  of  |  and  any  function  of  r\   are  also 
characteristic  coordinates.   It  is  not  necessary  to  find  £   and 
T)  explicitly  in  order  to  solve  the  equations  because  they  do  not 
appear  explicitly  in  the  equations.   However,  closed  form  solu- 
tions do  exist  and  in  Section  6  we  will  show  how  we  may  choose 
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characteristic  coordinates  i   and  r|  to  obtain  singular  solutions 
with  appropriate  infinities. 

Only  two  characteristics  i,    r\    pass  through  any  point  u,  v 
in  the  complex  hodograph  plane.   Thus,  we  often  denote  the  point 
by  |,  r) .   The  characteristic  curves  passing  through  a  point  on 
the  complex  sonic  line  have  the  same  characteristic  direction 
and  represent  a  singularity  of  the  characteristic  transformation. 
Therefore,  the  transformation  is  multiple-valued.   Two  charac- 
teristics £  and  T)  may  intersect  in  more  than  one  point  of  the 
hodograph  plane.   The  point  of  intersection  of  the  characteris- 
tics |  and  T)  will,  in  general,  depend  upon  the  paths  of  integra- 
tion of  the  equations  from  some  fixed  point.   In  Section  8  we 
explore  how  to  choose  paths  starting  in  the  subsonic  region  which 
give  rise  to  the  solution  in  the  supersonic  region. 
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5-  Incompressible  Flow  About  an  Airfoil 

We  assume  that  the  reader  is  already  familiar  with  many 
of  the  aspects  of  incompressible  flow  around  an  airfoil.   The 
purpose  of  this  section  is  to  formulate  the  problem  of  finding 
incompressible  flow  in  a  way  that  can  easily  be  generalized  to 
the  compressible  case.   In  particular,  we  find  a  characteristic 
coordinate  which  enables  us  to  express  the  singular  solution  as 
a  simple  pole.   Then  we  set  up  a  characteristic  initial  value 
problem  to  determine  a  regular  solution  which  can  be  added  to 
the  singular  solution  because  the  equations  of  motion  are  linear 
in  the  hodograph  plane.   Since  many  incompressible  solutions  are 
already  known,  it  is  possible  to  learn  the  relationship  between 
the  data  assigned  on  the  characteristics  and  the  resulting  air- 
foil.  This  knowledge  will  be  used  when  we  choose  data  on 
initial  characteristics  for  compressible  flow. 

We  start  by  applying  the  technique  of  the  previous  sec- 
tions.  Incompressible  flow  can  be  viewed  as  the  limiting  case 
of  compressible  flow  as  the  Mach  number  vanishes  or  the  sound 
speed  becomes  infinite.   In  either  case  the  equations  of  motion 
are  the  familiar  Cauchy-Riemann  equations 


(5-la)  ux  +  vy  =  ° 


(5-lb)  vx  +  uy  =  0 

Next  we  perform  the  hodograph  transformation  outlined  in  Section 
2  to  obtain 
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(5.2a)  yy+xu  =  0 


(5.2b)  -yu  +  xv =  ° 


which  we  note  is  also  the  limiting  case  of  (2.8a-b)  when  the 
sound  speed  becomes  infinite.   We  now  extend  u  and  v  into  the 
complex  domain  as  discussed  in  Section  ;>♦   We  compute  the 
characteristic  directions 

(5-3)  \  =   ±i 

which  can  also  be  found  by  letting  c  — >■  co  in  {J>.k). 

We  introduce  characteristic  coordinates  as  in  Section  4, 
and  we  find  that  equations  (4.7a-b)  and  (4.8a-b)  also  hold, 
where  A,  is  the  quantity  described  in  (5«3)«   For  the  readers 
convenience  we  rewrite  these  equations 

(5-^a)        u^  =  A+v^    ;    u^  =  \_v^ 

(5.4b)        y^  +A_x^  -  0  ;    y^  +^_^n   =  0 

We  turn  now  to  the  classical  theory  of  incompressible  flow 
about  a  body  with  unit  velocity  at  infinity  and  circulation  T. 
Any  such  flow  can  be  represented  by  a  complex  potential  x(z) 
z  =  x  +  iy,  where 


•V     1  V  /\ 

(5-5a)  x(z)  =  z+^  +  i  2?  loS  z 

z 


(5.5b)  z  =  H(z) 
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The  stream  function  if/   is  the  imaginary  part  of  x«   Equation 
(5- 5a)  represents  flow  with  circulation  r  about  a  unit  circular 
cylinder  in  the  auxiliary  z  plane,  while  (5* 5b)  is  a  mapping  of 
the  exterior  of  the  unit  circle  in  the  z  plane  onto  the  exterior 
of  the  body  defined  as  tj/  -   0.   The  complex  velocity  w  =  u  -  iv 
can  be  found  by  taking  the  derivative  of  x  with  respect  to  z. 
We  find  that 

(5.6)  w  =  5 i-  ;    K  =  J= 

H'(S)  47r 

First  we  confine  ourselves  to  the  case  of  a  circular 

<\ 

cylinder,    i.e.,    z   =   z,    so   that 


(5.7)  x(z)    =   z+^+21K   log   z 


and   thus 


,P   q,  ,        1    ,  2iK 

(5-8)  w  =  1  -  -£+-^- 


Equation    (5.8)    can  be   inverted   explicitly  and   yields 

(5.9)  z(w)  =  [(1  -  w-K2)l/2  +1K]"1 

2 

from  which  we  conclude  that  there  is  a  branch  point  at  w  =  1-  K 

and  a  pole  on  one  branch  of  the  w  plane  at  w  =  1.   In  the  general 
case,  we  see  from  (5.6)  that  the  inverse  function,  z  =  z(w),  can 
be  extremely  complicated  even  outside  the  body,  ip  >_   0,  z  ^  1. 
However,  for  a  wide  range  of  airfoils,  the  inverse  function  will 
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have  at  most  one  branch  point  of  first  order  outside  the  body, 
and  a  simple  pole  on  one  branch  at  w  =  1.   Let  us  introduce  the 
new  variable  r\ ,    defined  as 

(5.10)  ti(w)  =  (l-w-b2)1/2+  ib 

where  b  =  0  if  K  =  0  but  is  otherwise  arbitrary.   For  the  circu- 
lar cylinder  we  have 

(5-11)  z  =  - 

where  we  have  chosen  b  =  K. 

Consider  now  the  class  of  airfoils  which  can  be  expressed 
as 

(5.12)  z  =  i  +  fh) 

where  f  is  an  arbitrary  analytic  function  of  r\    outside  the  body. 

Then,  as  a  function  of  w,  z  has  a  first  order  branch  point  at 

2 

w  =  1  -  b   and  has  a  pole  on  one  sheet  at  w  =  1  for  all  solutions 

in  this  class.   Included  in  this  class  is  the  family  of  symmetric 
non-circulatory  elliptical  cylinders  defined  by 

(5.13a)  x  =  z  +  - 


(5.13b)  z  =  P(z+f) 

z 


where  E  is  an  ellipse  parameter  ranging  between  0  and  1,    and 


P  =  /1-E.   For  E  =  0  the  ellipse  is  a  circle  and  for  E  =  1  it  is 
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a  straight  line  segment.   Taking  the  derivative  of  x  with 
respect  to  z  we  find  that 

(5-14)  S  =  ^  (Et^  +  32)1/2 

1/2 
where  r\   =    (1-w)    '     ,    because   K  =   0  and  hence 

(5.15)  Z    =  4-    [ET!2+  p2)1/2   +  Et^  . 

Clearly  (5-15)  can  be  written  in  the  form  (5.12)  outside  the 
body,  since  from  (5.1*0  we  see  that  except  for  the  branch  point 

which  occurs  inside  the  body  at  z  =  0,  z{r\) is  regular. 

It  is  easy  to  see  that  w  is  a  characteristic  coordinate. 
It  then  follows  from  (5»10)  that  r\    is  a  characteristic  coordinate 
since  any  function  of  a  characteristic  coordinate  is  itself  a 
characteristic  coordinate.   Similarly,  we  can  show  that  %, 
defined  as 

(5.16)  i   =    (1-  w*-  b2)1/2-  ib  ;         w*  =  u  +  iv 

is  the  other  characteristic  coordinate.   From  this  we  conclude 
that  for  incompressible  flow  there  exists  a  characteristic  co- 
ordinate r\    such  that  the  solution  to  the  hodograph  equations  for 
flow  around  a  large  class  of  profiles  can  be  represented  as  a 
simple  pole  plus  an  analytic  function  of  t\   as  in  (5.12).   By 
allowing  f  to  have  branch  points,  we  may  further  extend  this 
class.   The  function  z  =  —  is  itself  a  solution  of  (5.4a-b)  and 
hence  the  class  of  solutions  which  represent  flow  around  an 
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airfoil  can  be  found  by  adding  any  regular  solution  of  (5.4a-b) 
to  — .   Prom  each  solution  of  an  initial  value  problem  or  a 
characteristic  initial  value  problem  of  (5.4a-b),  we  can  con- 
struct a  solution  for  flow  around  an  airfoil  simply  by  adding  x 
to  Re  (— )  and  y  to  Im  (— ). 

As  an  example,  the  solution  for  flow  around  an  elliptical 
cylinder  can  be  found  by  solving  the  characteristic  value  problem 


(5.17a)         x  =  g(|,r)o)   on  tj  =  r\Q 


(5.17b)         x  =  g(£Q,Ti)   on  i   =  iQ 

where 

(5.18)    g(C^)  =  1  (Er^+B2)1/2  +  E%   ,   -  i 

and  adding  the  real  and  imaginary  parts  of  x  and  y,  the  solutions 
to  this  characteristic  initial  value  problem,  respectively  to 
the  real  and  imaginary  parts  of  — .   It  is  easily  verified  that 
x  =  g(4,r|),  y  =  -igdjT])  is  the  solution  to  the  above  problem  so 
that  in  this  way  we  arrive  at  a  flow  about  an  elliptical  cylinder 
as  expected.   We  will  use  this  function  g  later  on  when  we  assign 
data  along  the  characteristics  in  the  compressible  case.   This 
problem  can  be  solved  numerically  by  the  technique  in  Section  7 
and  the  body  can  be  identified  by  examining  the  stream  function 
defined  by  equation  (2.5)-   In  fact,  this  was  done  to  check  out 
the  computer  program. 

By  our  construction,  it  is  clear  from  (5*12)  that  z  is 
single-valued  near  infinity.   However,  the  stream  function  if/ 
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must  also  be  single-valued  in  the  physical  plane  near  z  =  co 
To  show  that  the  stream  function  is  single-valued  in  the 
neighborhood  of  z  =  oo  we  must  show  that 


(5-19)  f  =   Im  (X)  =(p    wdz  =  0 

in  the  neighborhood  of  w  =  1.   If  we  carry  out  the  integration 
along  a  path  w  =  1+  e  ^,  0  <  <f>  <_  2tt  we  find  that  f   is  single- 
valued  provided  F  is  real.   For  symmetric  non-circulatory  air- 
foils, the  single-valuedness  follows  easily  from  the  symmetry 
of  the  integral. 

While  this  section  does  not  deal  with  an  equation  of  mixed 
type,  it  shows  us  how  to  obtain  a  solution  with  the  correct 
singular  behavior  by  finding  a  regular  solution  to  the  differen- 
tial equations. 
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6.  The  Singular  Solution 

In  Section  5  we  saw  how  a  specific  choice  of  the  charac- 
teristic coordinates  led  us  to  a  particularly  simple  form  of 
the  singular  solution  needed  to  represent  incompressible  flow 
about  an  airfoil.   We  noted  that  w  =  u  -  iv  was  a  characteristic 
coordinate.   The  general  solution  could  be  obtained  by  finding 
regular  solutions  to  the  equations  of  motion  and  adding  them  to 
the  singular  solution.   For  compressible  flow  the  procedure  is 
similar.   First  we  find  the  correct  form  of  the  singularity  for 
flow  around  an  airfoil.   Since  the  equations  are  linear  in  the 
hodograph  plane,  we  may  add  a  regular  solution,  found  by  solving 
a  characteristic  initial  value  problem,  to  the  singular  solution. 
In  analogy  with  the  previous  section  we  require  that  the  leading 
term  of  the  singular  solution  be  a  pole  in  the  appropriate 
characteristic  coordinate  r\.      In  order  to  form  a  solution  to  the 
equations,  we  need  a  lower  order  singular  term  and  a  regular 
term.   Thus,  we  expect  the  general  solution  to  have  form 


Z, 

(6.1)  z  =  _£  +  z2  log  T!  +  Z^ 


where  Z.,  I  =  1,2,5  are  regular  at  £  =  r\   =  0.   We  will  investi- 
gate the  conditions  necessary  to  impose  on  Z.,  and  Zp  in  order 
that  Z-,   is  the  solution  to  a  non-singular  differential  equation, 
We  shall  see  that  Z,  need  only  depend  on  i.      For  flow  with 
circulation,  Zp  is  closely  related  to  the  fundamental  solution 
and,  in  fact,  satisfies  equations  similar  to  those  of  the 
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Riemann  function.   For  flows  without  circulation,  Zp  can  be 
neglected. 

As  we  stated  in  Section  5*  incompressible  flow  could  be 
considered  as  the  limit  as  c  — ►  co  of  compressible  flow.   Thus, 
there  exists  a  characteristic  coordinate  s  with  the  property 
that  s  — ►  w  as  c  — ►  co  .   We  shall  endeavor  to  find  s  and  the  corre- 
sponding r)  in  this  section.   We  introduce  the  polar  coordinate  9 
in  the  hodograph  plane,  defined  as 

(6.2)  6   =  tan"1  - 

v    '  u 

A  solution  of  (4.6a-b)  can  be  found  in  terms  of  q  and  9   as 
illustrated  by  Sears  [21].   In  particular,  he  shows  that 

(6.3a)       f(q)  -9   -   const.   on  £   =  const. 

(6.3b)       f(q)  +9   =  const.   on  r\   =  const. 

where 

2  c2 

(6.3c)        f(q)    =  JL  sin-1[(7-D   \-7]  +|  sin^tfr+D  -|  -  7] 

^  c*  q 

Since  a  function  of  a  characteristic  is  itself  a  characteristic, 
it  is  clear  from  (6.3-b)  that  for  any  functions  F  and  G 

(6.4a)       s  =  F[f(q)  +  9]    -   const.   on  r\   =   const. 

(6.4b)        t  =  G[f(q)  -  9]    =   const.   on   ?  =  const. 
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We  want  to  choose  functions  F  and  G  in  such  a  way  that  s  and  t 
will  be  analogous  to  w  and  w*  in  the  incompressible  case.   We 
require  that 


(6.5a) 


lira  P[f (q)  +  0]    =   w  =  qe 
c— *■(© 


-i< 


(6.5b) 


lira  G[f (q)  -  0]    =   w*  =  qe 

C-*OD 


10 


We  will  show  that  the  function 


(6.6) 


F(r)  =  G(r)  =  Ce 


-lr 


has  the  desired  properties  for  s  and  t  for  an  appropriate  con- 
stant C.   With  this  choice  of  F  and  G,  s  and  t  assume  the  form 


(6.7) 


s  =  Ce-lf^e-i9  ;    t  =  Ce-if^ei0 


The  function  f,  expressed  in  terms  of  logarithms,  can  be 
written  in  the  more  convenient  form 


(6.8) 


f(q)    =   i    log   h(q) 


where 


(6.9a)      h(q)    =   Cq 


y/w. 


/ir^n/i7^ 


1-1/2 


(6.9b) 


II 


7 


7-1      2  ,_...,  x.2      _2. 


2 
c* 


q     ;      iP  =  (7+D ci  -qfy 


Writing  s  and  t  in  terms  of  the  function  h,  we  have 
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(6.10) 


s  =  h(q)e  1   ;    t  =  h(q)e1 


Finally,  by  choosing  the  constant  C  so  that  h(l)  =  1,  we  have 
the  desired  result. 

Continuing  our  analogy  with  the  previous  section,  we 
define  the  characteristics 


(6.11a) 

(6.11b) 


,-  ^2x1/2      .. 

•q=(l-s-b)/      +ib 


i  =  (i- 1  -bV/2  -ib 


and  seek  a  solution  to  our  problem  of  the  form  (6.1).   In 
particular,  we  assume  that  the  solution  can  be  written  in  the 
form 


(6.12a)   x  =  Re 


X  (|,il)  +X2(|,11)  log  n  +x^(i,r]) 


(6.12b)    y  =  Re 


Y  (^T))  +  Y2(^)  log  ti+Y5(|,t]) 


where  X1,    Y1,    i   =   1,2,3,    are   regular  at   i   =  r\    =  0.      Prom  the 
equations   of   motion    (4.8a-b)    we   find   that 


-3 


Y\  +  A_x| 


(6.13a)      Y^   +   A_X^  =  --^ *  -    (Y^+A_X^)    log  r) 


(6.13b)      Y^+A+X^  = 


,  ,         Y1  +  AX1         Y1  +  AX1  -  Y2  -  A   X2 

—£—  -  J ±3 fL  _  (Y2  +  A  Y2)log  t, 


In  order  to  be  able  to  compute  regular  solutions  we  require  that 
the  right  hand  side  of  these  equations  be  regular.   In  order  to 
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meet  these  requirements  the  following  conditions  must  hold  on 
the  characteristic  tj  =  0 

(6.14a)         Y*+A_xi  =  0  ;    Y1  +  A+XX  =  0 

(6.14b)         Y?+A  X?  =  0  ;    Y2  +  A,X2  -  0 

%  -  %  T\  +    T) 

(6.14c)  Y2  + A  X2  =  -(\J  X1 

2       2 

In  addition  to  these  conditions,  we  require  that  X  and  Y  be 

real  when  q  and  9   are  real  so  that  the  solution  is  single-valued 
in  the  real  physical  plane.   Finally,  we  set  X  (0,0)  =  1  so  that 
the  solution  corresponds  to  the  solution  for  the  case  of  in- 
compressible flow. 

Conditions  (6.l4a)  along  with  the  requirement  X  (0,0)  =  1 
yield  a  linear,  first  order,  ordinary  differential  equation  for 
X  whose  solution  can  be  written  in  closed  form 

(6.15)   XX(e,0)  =  expj   f    a(r,0)dr|;    a(£,ri)  =  -  A  _+^ 

Any  X   satisfying  (6.15)  and  the  corresponding  Y  can  be  used 

in  constructing  the  solution.   Thus,  we  may  choose  X  as  a 

2       2 

function  of  |  only.   From  (6.l4b)  we  see  that  X  and  Y  must 

satisfy  the  equations  of  motion  (4.8a-b)  plus  the  condition 
(6.l4c)  on  r\   =  0,  along  with  the  requirement  that  they  be  real 

at  real  points  of  the  hodograph.   This  means  that  we  can  find 

2      2 
X  and  Y  by  solving  a  characteristic  initial  value  problem  for 

the  equations  of  motion,  where  the  characteristic  initial  values 
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satisfy  a  linear  ordinary  differential  equation  on  t\   =  0.   We 

2  2 

see  that  X  and  Y  are  essentially  the  fundamental  solution  of 

112       2 
the  equations.   Having  found  X  ,  Y  ,  X  ,  and  Y  ,  the  corre- 

3  "3 

sponding  X^   and  Y^  can  be  found  by  solving  either  an  initial 

value  problem  or  a  characteristic  initial  value  problem.   In 
either  case  computational  difficulties  will  arise  in  the 
neighborhood  of  r\    =  0.   We  avoid  this  difficulty  by  assigning 
data  for  the  characteristic  initial  value  problem  along  charac- 
teristics 4  ,  T)   away  from  r\    =   0. 

The  stream  function  ip   can  easily  be  found  once  x  and  y 
are  known  by  integrating  (2.5).   We  find  that 


(6.l6)     i}/   =  —   /  (pvdx  -  pudy)  ;    ty  =   0  at   q  =  0 
P0  J 


In  this  way  the  body  is  defined  by  the  stream  line  f  =   0. 

In  order  for  the  solution  to  represent  flow  around  a  closed 
body,  it  is  necessary  that  if/   return  to  the  same  value  as  we  inte- 
grate on  a  closed  path  enclosing  the  body.   For  the  symmetric, 
non-circulatory  airfoil,  the  integral  around  a  closed  path 
vanishes  by  symmetry  arguments.   In  a  future  paper  on  non- 
symmetric  circulatory  flows,  a  condition  analogous  to  Re  r  =  0, 
for  incompressible  flow,  will  be  worked  out. 

For  the  symmetric  non-circulatory  case  a  further  simplifica- 
tion results  from  the  fact  that  s  is  an  even  function  of  r\.      From 
(6.11a)  we  have 


(6.17)  (V  )   =  (\J  ^  =   0  on  n  =  0 

+  T)  +  S  dT]  ' 
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p       p 

so  that  X  =  Y  =0  satisfies  (6.l4b-c).   Thus,  a  solution  to 
our  problem  can  be  found  simply  by  solving  a  characteristic 
initial  value  problem  for  the  non-homogeneous  system  of  equa- 
tions (6.15a-b)  with  X1  defined  in  (6.15)  and  X2  =  Y2  =  0. 
Since  the  choice  of  characteristic  initial  data  is  arbitrary, 
a  wide  range  of  airfoils  can  be  found  by  this  procedure. 
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7.  Finite  Difference  Method 

As  we  discovered  in  Section  6,  we  can  arrive  at  a  solu- 
tion of  the  equations  of  motion  for  flow  around  an  airfoil  by 

1   1 
solving  one  linear  ordinary  differential  equation  for  X  ,  Y  , 

and  one  system  of  non-homogeneous  linear  partial  differential 

"3   "3 

equations  for  X  ,  Y  .   Since  the  solution  to  the  ordinary 

differential  equation  is  known  in  closed  form  (6. 15),  the  first 
part  of  the  numerical  procedure  can  be  effected  by  almost  any 
quadrature  formula.   In  particular,  Simpson's  rule  is  a  simple 
and  rapid  procedure  to  accomplish  this,  although  care  must  be 
taken  to  avoid  the  region  around  which  A   is  equal  to  A_  as  the 
integrand  becomes  singular  there.   Instead,  we  choose  paths 
which  avoid  this  difficulty  as  explained  more  fully  in  the  next 
section. 

Having  determined  X  and  Y  ,  we  proceed  to  solve  a  charac- 
teristic initial  value  problem  for  the  equations 

(7.ia)    y|+a_x|  =  g1^^)  ;   g1(€,n)=  -  % 

Y  4-  i\   x 
(7.1b)     Y^+A+X^  =  g2U,ri)  ;    g2U,r))  =  ^— 

for  X^  and  Y^  derived  from  (6.15a-b)  for  the  symmetric  case 
without  circulation.   A  very  simple  method  of  integrating  these 
equations  has  been  developed  by  Massau  and  is  described  by 
Forsythe  and  Wasow  [5].   As  we  shall  see,  the  method  yields 
second  order  accuracy  and  will  always  work  away  from  the  sonic 
line. 
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Suppose  P-j^I-^t^)  and  P^^,^)  are  tw0  neart>y  points  in 
the  hodograph  plane  at  which  the  solution  is  known.   The  point 
P  (42,rj  )  in  the  hodograph  plane  where  the  characteristics  42 
and  ru  intersect,  as  well  as  the  solution  at  this  point,  can  be 
found  to  second  order  in  5,  where  6  is  the  distance  between  P, 
and  P„.   To  do  this  we  first  write  the  difference  approximations 

(7.2a)      u(P  )-u(P1)  =  A+(P1)[v(P3)  -  v(P1)] 

(7.2b)      u(P5)  -u(P2)  =  A_(P2)[v(P-5)  -v(P2)] 

for  equations  {H-.Ya.-b),    which  we  note  are  two  linear  algebraic 
equations  for  the  two  unknown  u(P-,)  and  v(P_).   However,  since 
A  and  A_  are  known  only  at  the  end  points  P.,  and  Pp,  it  is  easy 
to  see  that  u(P^)  and  v(P^.)  are  only  first  order  accurate  in  5. 
By  using  the  values  of  u(P^,)  and  v(P  )  obtained  in  this  way,  we 
can  find  A  (P-,)  and  A_(P^,)  to  first  order.   We  may  use  these 

values  to  compute  A  =  -p  [ A  (P  )  +  A  (P..  )  ]  and 

^    1 

A  =  ■£  [A_(P  ) +  A_(P_) ]  and  insert  them  in  (7.2a-b)  respectively 

to  obtain  u(P^)  and  v(P.,)  to  second  order  accuracy.   Having  done 

this  we  may  compute  g  =  2  [g  (P-,)  +  g  (P1 )  ] , 

g2  =  \   [g2(P3)  +g2(P2)l,  ^_  =  \   [A_(P3)  +A_(P1)],  and 

A+  =  \   [^+(p3)  +A_(P2)].   We  now  find  Y?   and  Y^  to  second  order 

in  5  by  solving  the  two  simultaneous  linear  algebraic  equations 

(7.3a)   Y3(P5)  -Y3(P1)  +A_[X5(P3)  -X3(P1)]  =  (^  -^Jg1   ' 

(7.3b)    Y5(P3)  -Y^(P2)  +  A+[X5(P3)-X2(P2)]  =  (r^-r^g2 
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which  result  by  applying  finite  differences  to  (7.1a-b). 

It  is  important  to  note  that  these  systems  of  equations 
involve  complex  quantities  and  must  be  solved  in  the  complex 
domain.   Thus  we  must  use  complex  arithmetic.   The  equations 
can  always  be  solved  provided  that  they  are  non-coincident. 
By  inspection  we  see  that  they  are  coincident  only  if  A  =  A  , 
i.e.,  we  are  on  the  complex  sonic  line.   If  we  are  near  the 
complex  sonic  line,  then  A   is  nearly  equal  to  A  and  the  equa- 
tions become  ill-conditioned.   Thus  we  can  expect  less  accuracy 
in  the  vicinity  of  the  sonic  line. 

For  the  characteristic  initial  value  problem  analytic  data 
are  prescribed  on  the  two  characteristic  surfaces  £,  =  £,  and 
r\   =  T]    .   On  the  two  dimensional  surface  £,  =  £,   we  pick  a  one 
dimensional  path  T)  =  t|(t)  starting  from  r\   =  t](0)  =  t\      and  put 
down  a  grid  r\    ,t)    ,t\    , . . .  ,r\    .      The  values  of  u  and  v  at  these 
points  can  be  found  explicitly  since  the  transformation  from 
characteristic  coordinates  to  the  hodograph  variables  is  known 
in  closed  form.   The  function  X^  is  prescribed  on  this  charac- 
teristic surface  and  Y^  is  found  integrating  equation  (7«lb). 
Similarly  we  lay  off  a  grid  on  a  path  i   =   £{o)    on  the  charac- 
teristic rj  =  rj   eminating  from  4  =  4   and  proceed  as  before  in 
obtaining  u,  v,  X  ,  and  Y  ,  except  that  equation  (7- la)  is  used 
in  place  of  [J. lb).      The  arc  length  parameters  o   and  t  are  real 
so  that  the  region  of  solution  can  easily  be  visualized  in  the 
( a, t) -plane.   We  see  that  the  domain  of  the  solution  is  the 
rectangle  bounded  by  the  four  characteristics  4  ,  r\    ,  i    ,  and  r\    . 
It  is  clear  that  the  method  of  Massau  can  be  applied  to 
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successive  points  in  each  row  starting  from  the  first  row  and 
working  up.   (See  Pig.  1) 

We  note  here  that  the  difference  equations  (7.2a-b)  for 
the  determination  of  u  and  v  are  used  only  to  facilitate  compu- 
tation.  It  is  of  course  possible  to  find  u  and  v  directly  from 
the  explicit  transformation  constructed  in  the  previous  section. 
In  solving  the  difference  equations  (7-2a-b),  A  and  A  are 
computed  at  the  midpoints  and  stored  for  later  use  in  the 
difference  equations  (7.j5a-b).   However,  along  the  initial 
characteristics  it  is  necessary  to  find  u  and  v  directly  from 
the  transformation. 

Physically  we  are  only  interested  in  the  solution  for 
real  values  of  u  and  v,  so  that  we  would  like  to  choose  paths 
on  the  initial  characteristic  surface  which  yield  as  many  real 
points  of  the  solution  as  possible.   By  examining  the  equations 
of  the  transformation  to  hodograph  variables,  we  find  that  if 
s  and  t  are  complex  conjugates,  the  point  u(s,t),  v(s,t)  is  in 
the  real  hodograph  plane  provided 


(7.4)  st  =  |s|2  =  |t|2  <   h2(cj 


From  equations  (6.11a-b)  we  see  that  £,  and  r\   are  complex  conju- 
gates when  s  and  t  are.   Thus,  by  choosing  paths  |(o),  t)(t)  which 
are  complex  conjugates,  we  find  that  the  points  of  the  two 
dimensional  grid,  shown  in  Fig.  1,  where  £,  =  r\   will  be  in  the 
real  hodograph  plane  provided 


(7-5)  |l-e2|  <   h(cj 
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For  these  conjugate  initial  paths,  we  obtain  the  solution  along 
a  path  in  the  real  hodograph  plane  starting  at  £.  "<]      and 
remaining  in  the  subsonic  region.   As  we  stated  earlier,  the 
real  part  of  x  and  the  real  part  of  y  are  the  desired  solutions 
along  this  path.   By  choosing  different  paths  along  £,  =  £   and 
•q  =  T)   we  find  the  solution  along  different  paths  in  the  real 
hodograph  plane.   We  see  that  any  complex  conjugate  pair  of 
paths  can  be  continued  only  until  |l  -  £.   =  h(c^.),  since  A   =  A_ 
there.   However,  by  choosing  conjugate  paths  which  cover  the 
initial  characteristic  £  =  £.   and  r\   =  r\      we  can  find  the  com- 
plete subsonic  portion  of  the  solution.   The  choice  of  paths  to 
find  the  solution  in  the  supersonic  region  will  be  discussed  in 
the  next  section. 

The  stream  function  f   can  be  found  by  using  the  second 
order  difference  scheme 

(7.6a)    V(P2)  =  ^(Px)  +p(P*)v(Pj[x(P2)  -  x(P1)] 

-  p(P*)u(P*)[y(P2)  -y(P1)] 


P2  -P1 

(7.6b)  P*  =  2~^ 


This  integration  need  only  be  carried  out  along  the  path  in  the 
real  hodograph  plane  along  which  the  solution  is  found.   For 
non-circulatory,  symmetric  flows  it  is  clear  that  if  we  choose 
£   =  T)   real,  then  v  =  0  and  hence  tyd    ,r\    )  =  0.   In  this  way, 
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the  body  can  be  identified  by  determining  the  points  where  t// 
vanishes  along  a  path  by  interpolation.   The  sonic  line  can 
be  determined,  but  to  less  accuracy,  by  taking  conjugate  paths 
which  terminate  when  A   =  A  . 
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8.  The  Sonic  Locus 

In  the  previous  section  we  showed  how  the  subsonic  por- 
tion of  the  flow  could  be  found  by  our  method.   At  the  sonic 
line  A   is  equal  to  A_  and  our  method  breaks  down.   The  inte- 
grand of  the  ordinary  differential  equation  needed  for  the 
singular  solution  becomes  singular  and  the  linear  difference 
equations  approximating  the  partial  differential  equations  be- 
come coincident.   These  singularities,  however,  are  singulari- 
ties of  the  characteristic  transformation  rather  than  singulari- 
ties of  the  flow,  and  hence  we  can,  in  general,  continue  our 
solution  into  the  supersonic  region.   If  the  problem  were  purely 
supersonic  in  nature,  clearly  we  could  solve  it  in  the  real 
domain  as  a  characteristic  initial  value  problem  by  our  method. 
In  this  case  A   and  A_  are  real  so  that  only  real  arithmetic 
would  be  necessary. 

For  conjugate  paths  along  conjugate  initial  characteristics 
£   and  r|   we  showed  that  we  could  continue  our  solution  until 

(8.1)         |i-n2l  =  |i-  e2|  =  h(c») 

We  shall  call  the  set  of  points  £,  on  r\    =  r\      and  the  set  of  points 
r\    on  |  =  |   which  satisfy  relation  (8.1)  the  sonic  locus.   For 
convenience  we  refer  to  a  point  i   on  r\   =  r\      as  belonging  to  the 
upper  sonic  locus,  while  points  t\    on  |  =  £   are  said  to  belong 
to  the  lower  sonic  locus.   Note  that  for  each  £,  on  the  upper 
sonic  locus  there  exists  a  point  r\   on  the  lower  sonic  locus  such 
that  r)  =  4-   The  solution  found  by  choosing  any  two  paths 
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terminating  at  these  two  conjugate  points  will  always  end  at  a 
point  on  the  sonic  line.   For  the  detached  shock  problem  a  more 
detailed  discussion  of  the  sonic  locus  can  be  found  in  a  paper 
by  Swenson  [25].   Note  that  here  the  sonic  locus  can  be  found 
explicitly  without  knowledge  of  the  solution,  since  the  problem 
is  linear  in  the  hodograph  plane. 

Let  us  for  a  moment  consider  the  characteristic  in  the 
supersonic  zone  of  the  real  hodograph  plane,  i.e.  consider  q 
and  0   real,  q  >  q  >_  c^.   From  equation  (6.3c)  we  see  that  f(q) 
is  real  so  that  from  (6.7)  we  have 

(8.2)  |s|  =  |t|  =  C 

and  thus  any  point  in  the  supersonic  region  is  the  intersection 
of  two  characteristics  s  and  t  each  having  magnitude  C.   As  a 
corollary  to  this  result,  we  can  conclude  that  for  every  point 
in  the  supersonic  region  characterized  by  the  two  characteristics 
s  and  t  passing  through  the  point,  the  corresponding  i   and  r\ 
defined  by  (6.11a-b)  with  b  =  0,  belong  to  the  upper  and  lower 
locus  respectively. 

We  can  see  this  in  another  way  which  does  not  rely  on  our 
explicit  knowledge  of  the  characteristic  transformation.   Through 
each  supersonic  point  in  the  hodograph  plane  pass  two  charac- 
teristics which  can  be  traced  back  to  the  sonic  line.   From  the 
sonic  line  we  can  trace  further  back  on  one  characteristic  until 

it  intersects  the  £   characteristic,  and  we  can  trace  the  other 

o 

further  back  until  it  intersects  the  r\      characteristic.   On 
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the  other  hand,  the  point  of  intersection  with  the  |   charac- 

'  o 

teristic  must  belong  to  the  lower  sonic  locus,  while  the  point 
of  intersection  with  the  r\      characteristic  must  be  on  the  upper 
sonic  locus. 

We  consider  a  path  on  £,  =    £   terminating  at  ( £,  ,t]  ),  where 
r|,  belongs  to  the  lower  sonic  locus,  and  we  consider  a  path  on 
r\   =  r\      terminating  at  (^,11),  where  £.,  belongs  to  the  upper 
sonic  locus.   We  further  assume  that  for  all  £,  and  all  t\    on 
these  paths,  A  ^  A  at  each  point  (£.,r|).   With  these  conditions 
imposed,  the  solution  can  be  found  numerically  by  the  computa- 
tional procedure  described  in  the  previous  section.   The  final 
point  of  solution  will  occur  at  the  intersection  of  the  charac- 
teristics 4-,  and  r\~  .   The  corresponding  s  and  t  must  have  magni- 
tude C  and  may,  in  fact,  be  a  point  in  the  supersonic  zone. 
However,  we  cannot  necessarily  conclude  that  the  point  (£,,t]  ) 
is  in  the  supersonic  zone  since  the  characteristic  transforma- 
tion is  not  single-valued.   By  choosing  the  paths  connecting  £. 
£,  and  r)  ,  r)   properly,  it  is  indeed  possible  to  end  up  on  the 
branch  of  the  transformation  such  that  the  point  (l,,^..)  belongs 
to  the  real  supersonic  zone. 

Two  questions  of  importance  naturally  arise.   First  of  all, 
can  paths  always  be  found  connecting  £,  £,  and  r\    ,  t^   with  the 
required  properties.   A  justification  for  the  existence  of  such 
a  path  is  that  the  set  of  points  for  which  A   =  A_  is  a  two 
dimensional  surface,  while  the  solution  space  for  any  two  initial 
paths  is  also  two  dimensional.   From  the  topology  of  the  complex 
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four-space  we  expect  that  we  can  pass  the  two  dimensional  solu- 
tion space  around  the  complex  sonic  line.   In  fact,  paths, 
similar  to  those  used  by  Swenson  [25]  for  the  detached  shock 
problem,  work  here.   In  order  to  describe  such  a  path,  we  con- 
sider the  point  ^M  on  the  upper  sonic  locus  which  is  the  complex 
conjugate  of  ru  .   We  assume  that  £„  is  to  the  right  of  £,  and 
thus,  r|  ,  the  complex  conjugate  of  £,  on  the  lower  sonic  locus, 
is  to  the  left  of  ru  .   To  find  the  point  [i-.,r\^)    in  the  super- 
sonic region  we  choose  a  path  on  the  t\      characteristic  starting 

from  4   and  terminating  at  £u  .   The  point  £LT  must  remain  to  the 
o  °  1  N 

left  of  this  path.   On  the  characteristic  £  we  begin  at  r\      and 
terminate  at  ru  keeping  r\      to  the  left  of  this  path.   If  we 
extend  the  path  £  ,  £   by  adding  the  points  £,...,£,  also  on 
the  upper  sonic  locus,  then  the  points  ( £p,r|  ) ,  ( £  ,t)    ) , 
.  ..,(£m,tu)  will  also  be  in  the  supersonic  zone  on  the  charac- 
teristic ru  .   In  this  way,  we  determine  the  solution  along  a 
characteristic  ru  in  the  supersonic  zone  starting  at  (£,,tu)  and 
ending  at  the  sonic  line  at  (£.  ,ru).   If  we  now  go  back  and 
extend  r\    ,  ru  to  the  point  r\^   on  the  lower  sonic  locus  we  can 
repeat  the  process  and  find  the  solution  along  the  characteris- 
tic T)    .      By  continuing  in  this  fashion,  it  is  possible  to  obtain 
the  solution  in  a  region  of  the  supersonic  zone  bounded  by  the 
sonic  line  and  the  two  characteristics  £-,  and  ru  with  one  set  of 
initial  paths  of  data.   The  region  of  solution  is  the  same  as 
would  be  obtained  if  we  had  solved  a  real  initial  value  problem 
with  the  initial  data  assigned  along  the  sonic  line.   In  Figs. 
2.a-c  we  have  illustrated  the  types  of  paths  involved  and  the 
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corresponding  region  of  solution  in  the  hodograph  plane. 

Secondly,  we  might  ask  whether  or  not  the  solution  is 
independent  of  path.   This  question  can  be  answered  affirma- 
tively only  if  we  know  that  there  is  no  singularity  of  the  solu- 
tion within  the  region  of  integration.   Thus,  we  might  end  up 
on  different  branches  of  the  solution  by  taking  different  initial 
paths  when  there  is  a  branch  point  in  the  solution.   In  this  way 
the  global  solution  can  be  obtained  and  branch  points  of  the  flow 
can  be  discovered.   The  next  section  is  devoted  to  a  discussion 
of  the  various  singularities  obtainable  by  this  method. 

Continuation  of  the  solution  to  the  supersonic  region  can 
be  achieved  with  one  set  of  initial  paths.   Moreover,  the  solu- 
tion is  obtained  along  characteristics  so  that  the  characteris- 
tics are  easily  plotted.   The  value  of  the  stream  function  is 
not  known  a  priori  anywhere  inside  the  supersonic  region.   Hence 
we  must  either  integrate  ip   in  complex  space  starting  from  (  £,  ,r|  ) 
where  if/  =   0,    or  we  must  start  the  integration  from  the  sonic 
line,  the  value  there  being  known  by  integrating  along  a  subsonic 
path.   The  latter  approach  has  been  adopted  in  the  numerical 
work  done  for  this  paper,  but  yields  less  accuracy  since  x,  y,  u, 
v  and  hence  ip   are  known  to  less  accuracy  at  the  sonic  line. 
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9.  Singularities  of  the  Flow 

Our  method  enables  us  to  find  airfoils  which  have  smooth 
flows  around  them.   The  only  types  of  singularities  of  the  flow 
that  can  be  computed  are  limiting  lines  and  branch  points. 
Smooth  flows  do  not  exist  for  every  airfoil.   In  general,  simple 
waves  and  shocks  will  form.   In  this  section  we  explain  how  the 
occurrence  and  location  of  limiting  lines  effects  the  behavior 
of  the  flow.   We  will  show  how  to  construct  a  solution  with  a 
weak  shock  when  the  limiting  lines  appear  outside  the  body  in 
the  decelerating  portion  of  the  flow. 

We  first  consider  what  types  of  singularities  may  occur 
in  the  flow  by  our  method.   In  Section  2,  we  showed  that  solu- 
tions for  which  the  Jacobian  j  vanishes  cannot  be  found  by  this 
procedure.   Thus  we  are  concerned  here  with  singularities  for 
which 


(9.1)  J=xy-xy=0 

v  ~ '  x '  uJ  v   vJ  u 


In  the  subsonic  region  the  Jacobian  J  never  vanishes  unless 

x  =v  =  x  =  v  =0.   This  can  be  deduced  from  the  relation 
u   Jv    v   Ju 

(9.2)     (c2-  u2)x^+2uvxuxy+  (c2  -  v2)x^  =  -(c2-u2)J 

obtained  from  equation  (2.8a),  which  shows  that  J  is  zero  only 

when  the  left  hand  side  of  the  equation  vanishes.   However,  for 

2    2     2     2 
subsonic  flow,  u  +  v  =  q  <   c  ,  and  the  left  hand  side  is  a 

positive  definite  quadratic  form  in  x  and  xv  so  that  the  result 

follows.   Moreover,  this  shows  that  the  Jacobian  can  vanish  only 
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at  isolated  points  in  the  physical  plane  and  hence  the  only 
singularities  possible  in  the  subsonic  portion  of  the  flow  are 
branch  points. 

For  supersonic  flow  J  may  change  sign  and  the  derivatives 
x  ,  x  -  y  ,  y  ,  need  not  all  vanish.   This  gives  rise  to  limiting 
lines,  which  are  discussed  in  some  detail  by  Courant  and 
Priedrichs  [h].      In  particular,  they  show  that  J  vanishes  along 
a  curve  in  the  hodograph  plane  whose  image  in  the  physical  plane 
is  a  fold  called  the  limiting  line.   They  show  further  that  this 
fold  is  an  envelope  of  one  of  the  families  of  characteristics, 
while  the  other  family  of  characteristics  cusps  along  it. 

It  becomes  evident  that  our  method  makes  it  easy  to  identi- 
fy limiting  lines  without  actually  computing  the  Jacobian  J. 
This  can  be  done  simply  by  plotting  the  characteristics  of  the 
flow  in  the  supersonic  region  of  the  physical  plane.   In  the 
previous  section  we  showed  how  the  characteristics  were  obtained 
as  a  biproduct  of  our  method  of  solution.   The  limiting  lines 
can  easily  be  recognized  by  examining  the  plot  because  one  family 
of  characteristics  cusps  there.   The  effect  of  a  limiting  line 
inside  the  body  can  be  shown  by  considering  the  speed  of  the  flow 
along  the  body.   As  the  limiting  line  approaches  the  body,  there 
is  a  more  rapid  variation  of  speed  near  the  limiting  line.   We 
illustrate  this  effect  by  plotting  the  pressure  coefficient  C 
defined  by 

p  "poo 
(9.3) 
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as  a  function  of  the  abscissa  x.   It  is  clear  from  equations 
(2.1)  and  (2.3)  that  C   is  a  monotonic  function  of  the  speed  q. 

IT 

Let  us  consider  the  case  where  the  limiting  line  is  inside 
the  body,  but  in  the  accelerating  portion  of  the  flow.   As  the 
limiting  line  approaches  the  body,  sharp  increases  or  peaks  are 
observed  in  the  pressure  distribution  curve.   A  region  of  rapid 
expansion,  similar  to  a  simple  wave,  is  formed.   This  type  of 
peaky  distribution  has  been  described  by  Pearcey  [20]  as  being 
a  good  design  for  shock-free  transonic  airfoil  sections.   For  a 
fixed  free  stream  Mach  number  the  limiting  line  can  be  made  to 
approach  the  body  by  increasing  the  thickness  to  chord  ratio. 
If  we  increase  this  ratio  too  much,  the  limiting  line  appears 
inside  the  flow,  making  it  impossible  to  physically  realize 
this  flow. 

The  effect  of  limiting  lines  in  the  decelerating  portion 
of  the  flow  is  reversed.   As  the  limiting  lines  approach  the 
body,  a  rapid  decrease  is  observed  in  the  pressure  distribution 
curve.   This  rapid  decrease  in  speed  would  probably  cause  the 
fluid  to  separate,  thus  changing  the  whole  nature  of  the  flow. 
We  would  expect  that  this  is  a  poor  aerodynamical  design  for  an 
airfoil.   One  advantage  to  such  a  solution,  however,  is  that  the 
supersonic  zone  occurs  further  back  along  the  wing  and  thus  a 
shock,  and  the  corresponding  shock-induced  separation,  would 
probably  be  delayed. 

When  limiting  lines  occur  in  the  decelerating  portion  of 
the  flow  outside  the  body,  it  may  be  possible  to  find  a  solution 
with  a  second  order  accurate  weak  shock  by  our  method.   To 
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illustrate  how  this  can  be  done,  we  consider  the  case  in  which 
there  are  two  limiting  lines  which  meet  on  the  sonic  line  out- 
side the  body  in  the  decelerating  portion  of  the  flow.   The 
solution  has  a  double  fold  in  the  physical  plane  and  can  be 
represented  on  three  sheets  joined  along  the  folds.   The  veloc- 
ity u,  v  as  a  function  of  x  and  y  will  be  three-valued  for  each 
point  x,  y  on  the  second  sheet  of  the  fold.   We  join  the  first 
sheet  of  the  solution  to  the  third  sheet  along  a  curve  on  which 
Prandtl's  relation 


(9.4)  q^  =   4 


is  satisfied,  where  q   is  the  speed  along  this  curve  on  the  first 
sheet  and  q,  is  the  corresponding  speed  on  the  third  sheet.   In 
this  way  the  jump  condition  for  the  normal  component  of  momentum 
is  approximately  satisfied  along  this  curve  (cf.  [4]).   In  order 
to  see  that  the  solution  obtained  in  this  way  represents  flow 
with  a  weak  shock  about  an  airfoil,  we  observe  that  the  stream 
function  t//   is  continuous  to  second  order  across  this  curve. 
Similarly,  we  can  ensure  conservation  of  tangential  momentum  if 
the  velocity  potential  <j),  defined  by 

(9.5)  <1>X  =  u  ;    <t>y  =  v 

is  continuous  in  the  physical  plane.   Indeed,  consider  the 
Jacobian  J  of  the  mapping  of  the  physical  plane  onto  the  (q>,V)- 
plane.   Prom  (2.5)  and  (9-5)  we  have 
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(9.6)        j  =  *xty-*y4>x  =  q2  £ 

and  thus,  to  first  order,  the  mapping  is  one-to-one  away  from 
stagnation.   We  can  conclude  that  f   and  <j)  are  continuous  to 
first  order  across  this  curve.   Thus  if  the  fold  is  not  very 
large,  this  curve  represents  a  shock  to  second  order  in  shock 
strength. 

In  this  way  we  can  generate  solutions  with  weak  shocks. 
We  can  view  the  appearance  of  a  weak  shock  as  a  perturbation  of 
smooth  transonic  flow  with  the  limiting  line  just  inside  the 
body.   To  explore  this  possibility,  a  flow  in  which  the  limiting 
lines  are  just  outside  the  body  in  the  region  of  deceleration 
of  the  flow  is  presented  in  the  next  section. 
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10.  Numerical  Results 

The  numerical  procedure  discussed  in  Section  7  has  been 
programmed  for  the  CDC  6600  at  the  AEC  Computing  Center  at  New 
York  University.   In  this  section  convergence  of  the  solution, 
the  time  and  storage  space  used  for  computation,  and  some  of  the 
resulting  airflows  are  discussed. 

As  we  pointed  out  in  Section  7,  the  finite  difference 
scheme  should  be  second  order  accurate  in  the  mesh  size  5.   To 
check  this,  we  ran  the  program  three  times  using  the  same  initial 
conditions  with  mesh  sizes  5  ,  5  /2,  and  6  /4  respectively. 
Table  1  shows  the  three  solutions  at  the  intersection  of  two 
given  characteristics  in  the  subsonic  portion  of  the  real  hodo- 
graph  plane.   The  fact  that  the  convergence  is  second  order  can 
be  seen  easily  by  examining  the  differences  between  the  three 
solutions.   The  final  column  shows  the  values  of  u,  v,  x,  y  and 
ip   by  extrapolating  these  three  solutions  to  zero  mesh  size.   We 
see  that  by  extrapolation  we  are  able  to  find  the  solution  to 
eight  significant  places.   Table  2  shows  similar  convergence  in 
the  supersonic  region.   Note  that  here  we  have  used  a  coarser 
grid  and  achieve  only  six  or  seven  significant  digits.   To 
verify  that  the  solution  is  consistent  with  the  original  equa- 
tions of  motion  in  the  real  plane  we  computed  the  solution  at 
five  points  (u  .v  ),  (u  ±h,v  ),  (u  ,v  ±h)  for  fixed  h  and 

r  v  O*  O  '   v  O     'o',xOO 

inserted  the  results  into  difference  equations  for  equations 
(2.8a-b).   The  results  of  this  check  are  presented  in  Table  3* 

The  program  was  written  so  as  to  accommodate  1000  mesh 
points  on  each  of  the  initial  characteristic  paths.   This 
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requires  about  40,000  sixty-bit  words  of  central  memory  because 
the  quantities  involved  are  complex.   A  typical  run  in  the  sub- 
sonic region  uses  200  grid  points  on  each  characteristic  and 
takes  about  20  seconds  to  compute  the  solution  accurate  to  five 
decimal  places,  provided  we  are  not  too  near  the  sonic  line. 
Near  the  sonic  line  the  solution  is  accurate  to  only  four 
significant  digits.   To  compute  the  whole  supersonic  zone,  a 
typical  initial  grid  uses  400  points  on  each  initial  charac- 
teristic and  takes  nearly  two  minutes  to  compute  the  solution 
accurate  to  four  or  five  places.   The  complete  solution  can  be 
found  and  plotted  in  about  six  minutes.   To  obtain  the  accuracy 
exhibited  in  Tables  1  and  2,  the  program  can  be  run  again  with 
a  different  mesh  spacing  and  a  deferred  approach  to  the  limit 
can  be  used. 

The  program,  as  written,  uses  an  unnecessarily  complicated 
form  of  the  singular  solution.   The  functions  X  and  Y  ,  in 
addition  to  satisfying  the  ordinary  differential  equation  (6.l4a), 
were  made  to  satisfy  the  original  differential  equations  of 
motion.   The  conditions  on  £,  =  0  are  taken  as  the  complex  conju- 
gates of  the  conditions  on  r\   -   0  found  by  solving  the  ordinary 
differential  equation.   Thus  it  is  easy  to  see  that  X  and  Y 
are  the  Riemann  functions  for  x  and  y.   The  initial  data  for  the 
regular  part  of  the  solution  is  assigned  along  the  characteris- 
tics I  =  £  and  n  =  ti  .   We  choose  £,  =  r\      real.   On  £   =  £     we 
^o      '    'o  o    'o  o 

choose  initial  conditions  of  the  form 

(10.1)        H1(t1)    =   g(£Q,T\)  -  eQ   log    (aQ-  t]  )  -  e1   log    [r\  -  a±) 
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where  g{i,r\)    is  defined  by  (5.18).   On  the  initial  characteris- 
tic r|  =  r]   we  take 

(10.2)   H2U)  =  g(^,Ti0)  -  eQ  log  (aQ-  t\q)  -  e1    log  [t)q  -  a±) 

In  this  way,  the  parameters  of  the  solution  are  M   ,  £,  ,  E,  e  , 

a  ,  e.,  and  an  ,  all  of  which  are  real.   Some  of  the  solutions 
o'   1      1' 

obtained  by  varying  these  parameters  are  plotted  in  Figs.  3-H- 
Since  all  of  the  flows  are  symmetric,  only  the  upper  half  of  the 
profile  is  plotted.   The  portion  below  the  axis  of  symmetry  is 
the  continuation  of  the  solution  inside  the  profile.   In  the 
supersonic  region  the  characteristics  are  plotted.   The  body  is 
shown  as  a  smooth  curve  connecting  over  100  points  of  the  solu- 
tion found  in  a  single  run.   The  sonic  line  can  be  found  by 
connecting  the  points  at  which  the  characteristics  cusp.   In  the 
subsonic  region  dashes  are  drawn  at  each  point  of  the  body  at 
which  the  solution  was  computed.   The  dashes  are  drawn  at  the 
angle  of  the  velocity  vector  through  that  point.   The  corre- 
sponding pressure  distribution  curve  was  computed  and  similarly 
plotted.   For  comparison,  all  of  the  airfoils  were  drawn  with 
the  same  chord  length. 

We  first  consider  the  case  where  £   =  e,  =  0.   From  Section 
5  we  see  that  as  M   — ►  0  the  solution  approaches  incompressible 
flow  around  an  elliptical  cylinder  with  ellipse  parameter  E, 
independent  of  %    .   Thus  we  expect  that  the  parameter  E  will  have 
a  similar  effect  for  compressible  flow.   As  the  free  stream  Mach 
number  M   is  increased,  the  body  becomes  thinner,  apparently 
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as  a  result  of  the  compression.   For  £   =  0  the  solution  has 

fore  and  aft  symmetry.   In  choosing  £   non-zero,  a  necessity 

for  computation,  the  airfoil  is  no  longer  fore  and  aft  symmetric, 

We  find  that  this  parameter  effects  where  the  limiting  lines 

appear.   For  £   positive,  a  limiting  line  occurs  near  the  back 

of  the  airfoil,  while  when  £   is  negative  the  limiting  line  is 

near  the  front.   For  £   =  0  there  are  limiting  lines  near  the 

front  and  back.   Figure  3  gives  an  example  of  flow  around  an 

airfoil  with  M   =  .7278,  E  =  .  381,  £Q  =  -A0.      The  airfoil  has 

a  thickness  to  chord  ratio  T/C  of  . 24l8,  and  the  maximum  Mach 

number  M    is  1.2560.   There  is  a  limiting  line  near  the  front 
max 

of  the  airfoil  inside  the  body  which  accounts  for  the  sharp  rise 

in  the  pressure  coefficient  in  that  region.   In  Fig.  4  we  have 

an  airfoil  with  M   =  .7906,  E  =  .  50,  £   =  -.32.   The  thickness 

00  0 

to  chord  ratio  for  this  profile  is  .1460,  and  1VI  v  =  1.2599. 
Similarly,  there  is  a  limiting  line  inside  the  body  in  the  for- 
ward region  of  the  flow,  and  the  resulting  pressure  distribution 
resembles  the  one  shown  in  Fig.  J>. 

Both  of  these  airfoils  would  probably  exhibit  poor  aero- 
dynamic properties.   The  speed  at  the  tail  falls  sharply,  and 
hence  separation  of  the  flow  from  the  airfoil  is  likely.   The 
separation  would  change  the  effective  shape  of  the  body  and  this 
might  lead  to  the  formation  of  shocks.   In  Fig.  ~},    the  speed 
returns  to  stagnation  even  though  there  is  a  corner  at  the 
trailing  edge.   The  angle  of  the  corner  is  60  because  the 
Jacobian  J  vanishes  at  the  stagnation  point  where,  as  we  see 
from  (9.6),  J  also  vanishes.   In  Fig.  4,  the  Jacobian  J  vanishes 
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away  from  the  stagnation  point  and  the  stream  line  cusps  at  an 

angle  of  0°.   The  coefficient  of  pressure  curve  does  not  have 

to  return  to  the  stagnation  value. 

In  order  to  improve  the  pressure  distribution  near  the 

trailing  edge,  we  insert  a  logarithmic  singularity  just  inside 

the  body  near  the  tail.   This  has  the  effect  of  stretching  the 

tail  so  that  the  speed  decreases  less  rapidly  and  moves  the  cusp 

so  that  the  speed  at  the  trailing  edge  is  greater.   The  pressure 

distribution  curve  is  less  steep  at  the  trailing  edge.   We  find 

that  this  tends  to  increase  the  curvature  at  the  nose,  and  in 

some  cases,  the  nose  becomes  cusped.   To  counter  this,  we  insert 

a  logarithmic  singularity  inside  the  body  near  the  nose  with  the 

opposite  sign.   The  parameters  e  and  e,  are  the  magnitudes  of 

the  logarithmic  singularities  at  the  tail  and  nose,  respectively. 

The  locations  of  the  respective  singularities  in  the  ^Q  plane 

are  given  by  the  parameters  a  and  a  .   For  the  flow  in  Fig.  5 

we  have  added  two  logarithmic  singularities  to  the  flow  in 

Fig.  4.   The  value  of  %     has  been  changed  to  bring  the  limiting 

line  closer  to  the  body  so  that  the  pressure  distribution  will 

be  peaky.   The  parameters  are  M^  =  .7906,  E  =  .50,  iQ  =   -A0, 

e     =    A,    a  =  .59,  e,  =  .15,  a,  =  -1.50.   Since  the  body  has 

been  elongated  near  the  tail,  T/C  decreases  to  .1279.   We  find 

M    =  1.2219.   It  is  clear  that  this  airfoil  would  be  less 
max        y 

likely  to  suffer  from  separation  effects  than  the  airfoil  in 
Fig.  4.   In  Fig.  6  we  have  an  airfoil  at  the  same  free  stream 
Mach  number  with  logarithmic  singularities  occurring  in  the  same 
place  in  the  hodograph  plane  as  in  Fig.  5.   The  magnitude  of  the 
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logarithm  term  at  the  nose  remains  unchanged,  while  the  magni- 
tude of  the  logarithm  term  at  the  tail  is  increased  to  .6.   The 
ellipse  parameter  E  has  been  decreased  to  .15  so  that  the  air- 
foil becomes  more  circular  (the  leading  edge  is  more  rounded) 
and  the  thickness  to  chord  ratio  is  increased.   The  initial 
characteristic  £      is  moved  to  -.16  to  keep  the  limiting  line  in- 
side the  body.   We  find  that  T/C  =  .1572  and  that  M    =  1.2708. 

max 

In  Figs.  7  and  8  we  have  chosen  £   positive  and  conse- 
quently the  limiting  line  appears  in  the  rearward  portion  of  the 
profile.   Both  flows  are  at  M   =  .7906  with  the  same  logarith- 
mic singularities  as  in  Fig.  5«   The  ellipse  parameter  E  is  .15 
for  both  airfoils.   In  Fig.  7  we  have  chosen  |   =  .14  and  the 
limiting  lines  are  just  inside  the  body,  while  £   =  .24  in  Fig. 
8  and  the  limiting  lines  meet  outside  the  body  in  the  manner 
described  in  the  previous  section.   We  explained  in  the  previous 
section  how  the  flow  in  Fig.  8  could  be  replaced  by  a  flow  with 
a  shock.   This  flow  can  be  viewed  as  a  perturbation  of  the  flow 
in  Fig.  7.   We  will  discuss  this  further  in  the  next  section. 

At  the  same  free  stream  Mach  number,  we  can  exhibit  a 
flow  in  which  there  are  limiting  lines  near  the  nose  and  the 
tail.   To  do  this,  we  choose  £,      close  to  zero  so  that  the  flow 
is  nearly  symmetric.   Figure  9  is  an  example  of  such  a  flow  with 
I  =  . 04  and  with  the  same  logarithmic  terms  as  used  in  Fig.  5- 
We  choose  E  =  0  to  bring  the  limiting  lines  close  to  the  body. 
This  makes  the  nose  very  blunt  and  increases  the  thickness  to 
chord  ratio  to  .1582.   We  find  that  the  speed  increases, 
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decreases,  increases,  and  decreases  again  as  we  go  from  nose 

to  tail.   The  maximum  Mach  number  is  1.2168. 

Finally,  we  present  two  flows  at  higher  free  stream  Mach 

numbers.   The  first,  in  Fig.  10,  is  at  M   =  .8392.   The  second, 

in  Fig.  11,  is  at  M   =  .8980.   We  use  the  same  logarithm  term 

that  has  been  used  in  Fig.  5  near  the  tail.   For  the  first  flow 

we  can  make  the  nose  rounded  by  inserting  a  logarithm  term  with 

en  =  .60,  a.,  =  -2.50.   For  the  flow  at  M   =  .8980,  we  were 
1      '   1  00 

unable  to  pick  e,  and  a   so  as  to  eliminate  the  cusp  at  the  nose. 
The  ellipse  parameter  for  the  first  flow  is  .40  and  the  initial 
characteristic,  |  ,  is  at  -.11.   We  find  that  T/C  =  .0956,  which 
is  thick  enough  to  be  practical.   For  the  second  airfoil 
T/C  =  .0465,  which  is  excessively  thin.   We  expect  that  the  upper 
limit  for  the  free  stream  Mach  number  needed  to  compute  realistic 
symmetric  airfoils  by  this  method  is  about  .86,  which  is  con- 
siderably higher  than  for  any  of  the  previous  methods. 
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11.  Physical  Interpretation  of  the  Results 

In  the  previous  sections  we  have  explored  a  method  to 
compute  airfoils  which,  In  theory,  have  shock-free  transonic 
flow  around  them.   The  question  arises  as  to  whether  these  air- 
foils would  exhibit  shock-free  flows  in  wind  tunnel  tests.   In 
the  past  twenty  years  a  controversy  has  arisen  as  to  whether 
shock-free  transonic  flows  could  exist  physically,  but  recent 
experimental  work  by  Koppe  [11],  Pearcey  [20],  and  Spee  and 
Uijlenhoet  [23]  indicates  that  essentially  it  does.   The  reader 
is  referred  to  papers  by  Nieuwland  [17]  and  Nieuwland  and  Spee 
[l8]  for  a  good  discussion  of  the  transonic  controversy. 

In  light  of  the  non-existence  theorem  of  Morawetz  [l6],  we 
see  that  a  profile  which  admits  a  shock-free  solution  is  an 
exception.   What  is  needed  is  a  well  posed  formulation  for  the 
boundary  value  problem  among  the  class  of  weak  solutions.   Weak 
solutions  satisfy  the  equations  of  motion  in  integral  form,  but 
allow  for  surfaces  of  discontinuity,  such  as  shocks.   Existence 
theorems  of  this  type  already  exist  for  certain  problems  in  fluid 
dynamics  (cf.  [19]).   With  such  a  theorem  we  would  be  able  to 
guarantee  that  the  solutions  found  by  our  method  are  meaningful, 
at  least  with  regard  to  the  mathematics.   A  perturbation  of  an 
airfoil  shape  found  by  this  technique  would  give  rise  to  a  solu- 
tion with  very  weak  shocks.   For  example,  the  smooth  airfoil  in 
Fig.  7  might  result  in  a  flow  with  a  weak  shock,  as  in  Fig.  8, 
if  the  body  is  slightly  distorted. 
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The  wind  tunnel  tests  done  in  Amsterdam  by  Spee  and 
Uijlenhoet  [23]  bring  out  many  interesting  facts  about  the 
stability  of  solutions  obtained  from  potential  theory.   First 
of  all,  the  solution  is  very  sensitive  to  the  shape  of  the  body. 
This  seems  to  account  for  the  difficulty  in  achieving  shock- 
free  transonic  flow  experimentally.   The  solution  must  be  calcu- 
lated very  accurately  and  the  airfoil  has  to  be  machined  with  a 
great  deal  of  precision.   We  have  shown  in  Section  10  the  capa- 
bility to  achieve  solutions  accurate  to  at  least  6  significant 
figures  by  our  method,  which  is  more  than  sufficient  for  con- 
struction of  a  virtually  shock-free  airfoil.   Spee  and 
Uijlenhoet  show  further  that  there  is  a  range  of  speeds  at  which 
the  airfoil  remains  nearly  shock-free  and  that,  if  the  free 
stream  Mach  number  is  approached  from  above  or  below,  shock-free 
flow  results.   Also,  for  a  small  range  of  the  angle  of  attack, 
the  solution  remains  almost  shock-free.   These  requirements  are, 
of  course,  necessary  if  we  are  to  build  aircraft  which  fly  in 
the  transonic  regime  without  shocks  and  the  resulting  wave  drag. 

Kuo  [13]  gave  arguments  to  show  that  shock-free  transonic 
flows  would  be  physically  unstable.   His  method  was  based  on  the 
amplification  of  accoustic  waves  in  transonic  potential  flow  near 
the  sonic  line.   Essentially,  Kuo ' s  argument  is  an  elaboration 
of  the  well-known  fact  that  shock  wave  has  no  equilibrium  posi- 
tion in  a  one  dimensional  decelerating  flow.   His  arguments  do 
not  take  the  two  dimensional  nature  of  this  problem  in  account. 
Spee  [22]  has  used  a  two  dimensional  argument  to  show  that  these 

flows  are  stable  provided  that  M    <  I.581. 

^  max  — 
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Without  experimentation  we  cannot  predict  what  kind  of 
success  to  expect  from  the  airfoils  presented  in  this  paper. 
In  light  of  the  experiments  in  Amsterdam  by  Spee  and  Uijlenhoet, 
we  expect  similar  success  with  our  airfoils  even  though  ours 
have  a  larger  supersonic  region.   We  believe  that  the  airfoils 
in  Figs.  3  and  4-  would  result  in  separation  of  the  boundary 
layer.   Airfoils  5>  6  and  10  appear  to  be  better  candidates  for 
physical  reality  in  view  of  the  experiments  of  Pearcey  [20]  in 
which  he  shows  that  a  peaky  distribution  in  the  accelerating 
portion  of  the  flow  is  a  good  design  for  shock-free  wings.   The 
flow  in  Fig.  7  might  also  prove  to  be  a  good  design.   The  super- 
sonic zone  occurs  closer  to  the  tail,  and  thus,  if  a  small  shock 
were  to  occur,  the  shock-induced  separation  would  appear  at  a 
point  further  back  along  the  wing.   This,  of  course,  is  specula- 
tion and  no  conclusions  can  be  drawn  without  wind  tunnel 
experiments. 

We  conclude  by  noting  that  we  have  constructed  a  so-called 
black  box  which  generates  airfoils  which  have  shock-free  tran- 
sonic flow  about  them  in  a  matter  of  some  four  or  five  minutes. 
For  a  given  thickness  to  chord  ratio  we  are  able  to  find  flows 
at  a  higher  free  stream  Mach  number  than  has  been  obtained  by 
other  methods.   This,  of  course,  means  that  planes  designed  by 
this  technique  would  fly  faster.   We  have  control  over  the  type 
of  pressure  distribution  which  results  so  that  we  are  able  to 
find  flows  which  are  good  candidates  for  physical  reality.   Our 
method  has  a  very  natural  extension  to  flows  with  circulation. 
This  will  be  the  topic  of  a  future  paper. 
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5 

.01 

.005 

.0025 

Extrapolation  to 
0 

u 

.8121013579 

.8121013774 

.8121013821 

.81210138 

V 

.4859075293 

.4859077488 

.4859078011 

.48590782 

X 

-2.6887276841 

-2.6887212866 

-2.6887196846 

-2.68871915 

y 

.0960896118 

.0960930475 

.0960939105 

.09609420 

i> 

-.0390955281 

-.0390990577 

-.0390999385 

-.03910023 

Table  1 


Convergence  at  a  subsonic  point 
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6 

.02 

.01 

.005 

Extrapolation  to 
0 

Re  u 

1.8225435832 

1.8226246466 

1.8226456682 

I.8226528 

Im  u 

-.0000477829 

-.0000127384 

-.OOOOO32881 

-.0000000 

Re  v 

.4248831668 

.4250125398 

.4250465989 

.4250583 

Im  v 

-.0000547324 

-.0000146533 

-.0000037767 

.0000000 

X 

1.0242850615 

I.0253285855 

1.0256005600 

1.025693 

y 

-I.4475380824 

-1.4471360304 

-I.4470267282 

-1.446988 

i, 

.1219670856 

.1222614743 

.1222958932 

.12230 

Table  2 


Convergence  at  a  supersonic  point 
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u 

5 

2   2 

c  -u 

=   .80 
=   .0025 
=  1.00 

v  =  .40 

q2  =  .80 

2uv  =  .64 

h  =   .0005 
c2  =  1.64 
c2-v2  =  1.48 

x ( u+h , v ) 
x(u-h,v) 
x ( u , v+h ) 
x(u, v-h) 

=  -2.7058917^23 
=  -2.7065876332 
=  -2.7061085400 
=  -2.7063723634 

y(u+h,v)  =  .1820240630 
y(u-h,v)  =  .1817602438 
y(u,v+h)  =  .1812935915 
y(u,v-h)  =  .1824923545 

X 

u 

•v  x(u+h,v)  -x(u-h,v) 
2h 

.6958909 

xv 

~  x(u,v+h)  -  x(u,v-h) 
2h 

.2638243 

yu 

~  y(u+h,v)  -y(u-h,v) 
2h 

.2638192 

yv 

~  y(u,v+h)  -y(u,v-h) 
2h 

-I.1987630 

-yu  +  xv  =  .0000051 

(c2. 

-u2)yv  +2uvxy+  (c2-v2)xi. 

=  .0000025 

Table  3 

A  check  to  show  consistency  of  the  solution  with 
the  original  differential  equations. 
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Pig.  1 

A  typical  grid  in  the  subsonic  portion  with  11  points  on 
each  initial  characteristic.   The  point  with  the  arrows 
pointing  towards  it  will  be  the  next  point  computed. 
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Fig.  2a 

Typical  initial  paths  along  the  initial  characteristics 
T}=T]  (top)  and  £  =  £  (bottom)  which  yield  the  solution 
in  a  region  of  the  supersonic  zone. 
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Fig.  2b 

The  corresponding  grid  for  the  paths  of  initial 
data  in  Fig.  2a  shown  in  the  o-\±   plane. 
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Fig.  2c 

The  region  of  solution  in  the  real  hodograph  plane 
found  by  using  the  initial  paths  of  Figs.  2a  and  2b 
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Flow  at  Mqo  =  .7278,  E  =  .581,  iQ   = 
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Flow  at  M^    =    .7906,    E   =    .50,    £ 
Wv  =   1.2399   and   T/C   -    .1460. 
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Fig.    6 
Flow  at  M^    =    .7906,    E   =    .15,    iQ   =    -.16,    e0    =    .6,    aQ   =    .59, 
en    =    .15,    an    =    -1.50.      M   o      =   1.2708   and   T/C   =    .1372. 
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Fig.    7 
Flow  at  M^    =    .7906,    E   =    .15,    CQ  =    .14,    eQ   =    .4,    aQ   =    .59, 
e     =    .15,    a     =   -1.50.      M  =   I.2038  and  T/c  =    .1^75- 

-1-  -L  UldiA. 
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Fig.    8 
Flow  at  IVI        =    .7906,    E   =    .15,    £      =    .2h,    e      =    A,    a     =    .59, 


00 


o 


o 


e1   =    .15,    a1   =   -1.50.      Mmax  =   1.2957  and   T/C   =    .1501. 
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Fig.  9 
7906,  E  :=  0,  ^  =  .04,  eo  =  A,    aQ  =  .59, 
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■1.50.   M    =  1.2168  and  T/C  =  .1582. 
^     max 
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Fig.    10 
.8392,    E   =    .4,    £      =    -.11,    e      =    A,    a     =    .59, 
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■2.50.      M  =    1.1978   and   T/C   =    .0956. 

max 
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1.1204   and  T/c   =    .0465. 
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FORMAT  (5F5.2) 

END 


SUBR 
THIS 
THE 
REGI 
COMP 
COMM 
COMM 
COMM 
COMM 
DATA 
1/(0, 
C(l> 
C(2> 
C(6) 
C(7> 
C(B) 
C(5> 
C(9) 
C(10 
C(ll 
C(l2 
C(16 
C(19 

CALL 
C(l3 
C(l4 

C(l5 
C(l7 

QST 
UO  ■ 
C  s 
C(4> 


OUTIN 
ROUT 
PROGR 
ON  FR 
LEX  U 
ON  /A 
ON  /9 
ON  /C 
ON  /D 
ARQi 

il,  )> 

=  QA 
s  GA 
=  -S 

8   C( 

»  C( 
s  C/ 
.5 
1 


2 

C 

c 

HSM 
■  C 


E 

■  s 

s  QST 

C(3> 

2,/(G 

=  RH 


E  CONS 

INE  CO 

AM  AND 

OM  W  a 

0,U»TA 

/  C(20 

/  U0(4 

/  ALP< 

/  N,M, 

ARG2(1 

3*d.# 

MMA-1, 

MMA+1, 

QRT(C( 

D/QST 

2)*QST 

C(8> 

*C(1) 

,/C(l> 

,25*C( 

,*C(3> 

(6)*C( 

(3)*C( 

(C(19) 

(19)/C 

5*(GAM 

5*C(9) 

QRT(1, 

*C(14) 


T  (H) 

MPUTES  AND  STORES  THE  CONSTANTS  WHICH  ARE  USED  IN 

FINDS  THE  SINGULAR  PART  CF  THE  SOLUTION  IN  THE 

1  TO  THE  INITIAL  CHAR ACTER  I  S I T I CS 
0,ALP,BET,B,H,B2,ARG,ARG2,V 
),GAMMA,QST,ARG(5),ARG2(5) 
,1024),V(4i1024),U(2,900)#TAO(2,10  24) 
1024),BET(1024),B,B2 
IN.MM.NRN 

),C(l3),C(18),C(20),UO(2),U0(3),BETd),ALPd)»UO(4) 
O,),2*(0.*l,)»l.,,5,0,f4*(0,#0.),(1..0.)/ 


2)/C) 


2)*QST 

7) 
3) 
iC(l3),C(l)> 

(13) 

MA*1.) 

♦C(14) 

/(C(14)*QST/C(3)-C(9) ) ) 


AMMA*RHO(UO) ) 
0(UO)*(QST-C(9)*C<3)*C(3)> 
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CALL  TOFH  ((0 

B2  =  TA0(2> 

B  =  2,*R2 

DO  10  J  =  2.N 

8ET<J>  =  BET(J-D*H 

ALP(J)  *  CONJR(BET(J)> 

CALL  INVAL  (2.N.BET) 

CALL  CONJ  (4*N,UO.V) 

CALL  XAri  (2.1. N.N.I) 

CALL  CONJ  (4.U0.V) 

NM  s  N-l 

DO  20  J  =  2.NM 


0. ), UO.TAO(l)  ) 


ALP( J) 
BET< J> 
RETURN 

END 


TAOCli J) 
TA0(2, J) 


C 
C 


10 


20 


30 


40 


50 
60 


SUBROUTINE  AANDB  (  N  .  ti'4,  ALP  .  X  ,  L  .  MPS  I  ,  MT  ) 

READS  IN  THE  PATHS  ON  THE  COMPLEX  CHARACTERISTICS  AND  SETS  UP  THE 

INITIAL  GRID 

COMPLEX  ALP(1),H 

LBM  =  1024 

READ  50.  M.MT 

MPSI  s  K  =  N 

DO  40  I  =  l.M 

READ  60.  ALP(LBM),LL 

MM  =  CARS(ALP(LBM)-ALP(»<)  )/X*,999 

MM  =  MM*LL*L 

KK  =  K*l 

IF  (MM)  20,10,30 

MM  s  1 

GO  TO  30 

MM  s  -MM 

MPSI  =  *>MM 

H  s  (ALP(LBM)-ALP(K) )/FLOAT(MM) 

K  =  K+MM 


DO  40  JJ 
ALP( JJ) 
NN  =  K 
RETURN 
FORMAT 
FORMAT 

END 


i  KK.K 
ALP(  JJ- 


1)*H 


(213) 

(2E20 


12,15) 
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SUBROUTINE  INVAL  (K.N.BET) 

COMPUTES  U  AND  V  ALONG  THE  CHARACTERISTIC  ETA 

COMMON  /B/  UO(4,1024,2),U(2,900),TAO<2,1024) 

COMPLEX  H,TA0,BETA,U,UC,BET(1),T<8) 

DO  10  I  =  K,N 

IMi  s  1-1 

H  s  BET( I  >-BET( IMi) 

T(7)  =  TAO(l, I ) 

=  BET< IMI) 

TOFH  (BETA*, ?5*H,U0(1, IMI), T) 

TOFH  (HETA*,5*H,U0<lj IM1),T(4) ) 

=  BETA+H 

TOFH  (HETA,UO(l, IM1),TA0<1, I ) ) 

XYO  <H/6,,TA0<1, IM1),T, I ) 


BETA 
CALL 
CALL 
BETA 
CALL 
CALL 
RETURN 

END 


C 
C 


SUBROU 
COMPUT 

charac 

COMMQN 
COMPLE 
IMI  * 

Y  =  SI 
X  =  UO 

Y  =  SI 
U0(4, I 

Y  z  SI 
U0(3,  I 
RETURN 
END 


TINE  XYO  <H,T,S, I > 

ES  X  AND  Y,  THE  COEFFICIENT  OF  THE  SINGULAR  TFRM,  ALCNG  THE 

TERISTIC  ETA  «  0, 

/B/    UO(4,1024,2),U(2,900),TAO<2,1024) 
X  H,T(l),S(l>,UO,U,TAO,X,Y»SIMP 

1-1 

MPCS(7),S<3),S(6),H/2, ) 

(4,  IM1)*CEXP(-Y)*S(4) 

MP(S(7),S(6),T<5)#H> 

)  =  UD(4, IM1)*CEXP(-V) 

MP(UO(4,lMl)*T*S(7),X*S(6)iU0(4,I)*T(3)*T(5),H) 

)  «  U0(3, IM1)*Y 


C 

C 
C 


SUBROUTINE  TOFH  <B,U,TT> 

GIVEN  THE  VALUE  OF  ETA,  THIS  ROUTINE  COMPUTES  LAMBDA*  ANC 

AS  WELL  AS  THE  FUNCTION  A  WHICH  IS  USED  IN  THE  OUADRATLRE 

ALONG  THE  CHARACTERISTIC  THROUGH  Hal, 

COMMON  /A/  C(20), GAMMA, QST,ARG(5) 

COMPLEX  B,U(6),TT(4)#T,CS,P(4),D(4),TQ,TTH,QS,HP,ARG,AR 

T  r  C(3)-B*R 

CALL  VELOC  (T,U,QS.HP) 

CS  =  QST*C(9)*QS 

P  s  CSQRT(CS*(QS-CS) ) 

CALL  BRAN  < P, ARG ( 5  )  , AR, K  ) 


L AMRDA- 
FORMULAS 
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IF  (K)  20, in, 20 

PM)  =  U(6)*U(6> 

ARG<5)  =  AR 

P(2)  =  U(5)*U<6> 

DM)  =  -U(5)*U(5) 

P(3)  ■  CS  +  DM) 

D<1>  =  2,*P(2> 

D(2>  =  -C(9)*D(4)/QS 

D(3>  =  QS«2.*PM> 

TT  =  (P<2)*P)/P(3) 

P<4)  =  (C(14)*CS*C(1D)/P 

TT(2)  s  (P(2)-P)/P(3) 

D<4)  =  P(2)/QS-P(4) 

TQ  =  (D(4)-TT(2)*D(2))/P(3) 

TTH  =  (D(3)-TT(2)#D)/P<3) 

TT<3)  =  <TQ*C<12)/HP*TTHMO,,l 

RETURN 

CALL  EXIT 

END 


)/T)*B/(TT-TT(2) ) 


C 

c 
c 


10 
20 
30 

40 


SUBROUTINE  VELOC  (B,U,CS,HP) 

TRANSFORMS  CHARACTERISTIC  VARIABLES  TO  HCDOQRAPH  VARIABLES 

GIVEN  T  AND  THE  VELOCITY  AT  THE  PRECEEDING  POINT,  THE  VELOCITY  IS 

COMPUTED  ALONG  THE  INITIAL  CHARACTER  I S I TC  S  «  If  AT  T  ■  T 

COMMON  /A/  C(20>, GAMMA, QST,ARG<5) 

COMPLEX  A,B.U(6),QS,H,HP#H1,W,WS,D,QG,ARG,AR 

DATA  ER/.1E-11/ 

D  =  b*C(3) 

QO  =  U(1)*U(1)*U(2)*U(2> 

DO  10  I  s  1,50 

CALL  HSH  (QO,H,HP) 

hi  =  H-n 

QS  =  UO-H1/HP 

IF  (CABS(OS-QO) .LE.EK)  30,10 

QO  =  QS 

PRINT  50,  I 

PRINT  60,  W,B,QS,Q0,ARG(3),AR 

CALL  EXIT 

W  =  CSQRT(QS*B/C(3) ) 

CALL  BRAN  ( W , ARG ( 3 ) , A R , K > 

IF  (K)  20,40,20 

ARG<3)  =  AR 

WS  =  C(3)*W/B 

U<5>  =  ,5*(W*WS) 

U(6)  =  (0., ,5)*(W-WS> 

RETURN 
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FORMAT 
FORMAT 
END 


(  110) 
(4E20 


11) 


SUBROUTINE  HSH  <Q,H,HP) 

GIVEN  U**2+V**2,  THIS  ROUTINE  COMPUTES  THE  FUNCTION  H, 

DERIVATIVE,  WHICH  IS  THE  PRODUCT  OF  THE  CHARACTERISTIC 

S  AND  T, 

COMMON  /A/  C(20>, GAMMA, QST, ARG(5) 

COMPLEX  Q,H,HP,TH1,TH2,S1,S2,T,AR,ARG 

THl  s  GAMMA-C(7)*Q 

TH2  -    C<8)-GAMMA*0 

51  =  CSQRT(TH1*TH1-1. ) 
CALL  BRAN  <S1,ARG(4),AR,K) 
IF  (K)  10,20,10 

10  PRINT  40,  C, GAMMA, QST,ER#ARG,AR,Q,H, HP 

CALL  EXIT 
20  ARG<4)  =  AR 

52  =  CSQRT(TH2*TH2-Q*Q) 
CALL  BRAN  ( S2 , ARG ( 2 ) , AR, K ) 
IF  (K)  10,30,10 

30  ARG(2)  *  AR 

T  s  CEXP(C(6)*CL0G<TH1-S1) ) 

T  s  C(13)*T/(TH2+S2> 

H  a  G)*T 

HP  =  C<16)*H/S1*C(8)*T/S2 

RETURN 
40  FORMAT  (10F12.7) 

END 


AND  IT*S 
COORDINATES 


COMPLEX  FUNCTION  S I  MP < A , 8 , C , H > 
INTEGRATION  BY  SlMPSON*S  RULE 
COMPLEX  A,B,C,H 
SIMP  B  H*(A*4.*B*C) 
RETURN 

END 


SUBROUTINE  STEP  (AA.NN) 

THIS  IS  THE  MAIN  CALLING  ROUTINE  FOR  FINDING  THE  SOLUTION 

COMPLEX  UO,U,TAO,ARG,ARG2,ALPiBET,V 


PAGE 


10 


20 


30 


COMMON  /A/  C(20>, GAMMA, QST, ARG<5) ,ARG2(5) 

COMMON/M/UO(4»in24)tV(4»1024),Li(2#900)»TAO<?»10  24),D(600)»A(4«600) 

COMMON  /C/  ALP(1024) ,H£T(1024) 

COMMON  /[)/     nj  ,  M ,  IN.MT.NfiN,  I  I  ,  JJ,KK 

DATA  KQ.LBM/0, 1(124/ 

N  P 1  s  N  *  1 

NMl  =  N-l 

ARG  s  ARG2<5) 

CALL  KAHii  (ARR2<2)»AR0(2>  ) 

D  =  0, 

M  s  0 

CALL  GAPG  ( V, UO  < 1 . N  )  ) 

CALL.  TOKH  (PET(N)  ,UO(l#N)  ,TAO<l»N)  ) 

ALP(N)  =  BET(N) 

CALL  AANUB  ( N » NN, ALP » A  A #  I  N ,  I  I , KK  ) 

CALL  INVAL  (NP1,NN,ALP) 

DO  10  I  =  2, NMl 

TAO<l»  I  )  =  ALP(  I  ) 

TA0<2.  I  )  =  BET<  I) 

NO  =  NN-N 

CALL    CONJ    (MQ,ALP(NPD#ALP<NP1)  ) 

CALL  CONJ  (4*NQfUO(lf NP1 ) , V < 1, NP1  )  ) 

CALL  CONJ  (4*NM1,V(1#2)#U0(1#2>) 

CALL  TAU  (UO(l,N),TAO(l#N>  ) 

CALL  XAB  (NP1,1,NN»N,1) 

ARG2  =  ARG2(5) 

IF  (KK)  20#20i40 

CALL  CONJ  (4*NQjV(l.NPl)#UO<l»NPD) 

CALL  CONJ  (4,V(l#N),UO(l#N) ) 

DO  30  I  =  N.NN 

CALL  TAU  (UO(l,  I).TAO(l,  I  )  ) 


40 


JJ  = 
MT  = 
CALL 
CALL 
CALL 
WRITE 


II 
NN 
GUO 
SAVE 
XAB 
(1) 


(NN+l-N.U, ALP(N) ) 
<UO(l,N),U0<3#N)#U,N,N) 

(NP1,N,NN,NM,2) 
K0#M, ( (A( If J)#  I  -    1»4)#D< J)i J 


i#m) 


GO  TO  60 

CALL  AANDB  (N, MM, RET, AA# IN, JJ.KK) 

CALL  GARU  < ARG2 ( 2 > , ARG < 2 ) > 

CALL  GAWG  (V.UO(l.N)) 

CALL  TOFH  (HET(N) ,UO(l#N),TAO(ljN) > 

CALL  IMVAL  (NP1,MM,BET) 

CALL  GUO  (MM*1-N,IJ,  ALP(N)  ) 

MT  s  MM 

CALL  XAB  (2,N, I  I, MM, 2) 

WRITE  (1)  KK,M, < <A( I, J>« J  =  1.4),D(J),J  - 


l.M) 


D(400)  = 
MTT  s  MM' 


D(M) 
JJ 
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DO  50  I  =  l.MTT 

M  s  0 

IP  s  I I+I 

MT  s  MM- I 

CALL  XAB  (  IP.N,  JP.MT.2) 

WHITE  (1)  KK.M, ( (A(K, J)»K 

CALL  HEAD 

RETURN 

END 


l»4),D(J),w  *  1,M) 


SUBROUTINE 
THIS  ROUTI 
TIAL  CHARA 
COMPLEX  UO 
COMMON  /A/ 
COMMON  /B/ 
COMMON  /C/ 
COMMON  /D/ 
READ  (2)  ( 
M  s  K 
NP1  ■  N*l 
NNPl  =  NN 
NQ  »  NN-N 
ALP(NN)  s 
CALL  CONJ 


STP2  <AA,NN,K) 
NE  ENABLES  US  TO  BRANCH  OFF  INTO  MANY  PATHS  CN  Tl-E  IM 
CTERISTIC 
,U,TAO,ARQ»ARG2,ALP.B(=T,V 

C(20)*QAMMA,(JST,ARG(5),ARG2(5) 

UO(4,1024)iV<4,1024)fU(2.900)i 

ALP(1024).RET(1024) 

N»M,  IN.N1.NRN,  I  I  , JJ,KK 
(UO<  I, J),  I  =  1*4), (TAO( I 


TAO(2,1024),D 


J)  ,l<  I. j)#  I 


1,2) 


If  NM) 


10 


*1 


CALL 
CALL 
CALL 
CALL 

N2  ■ 


CONJ 
TOFH 
AANDB 
INVAL 

Nl-NN 


ALP(NN)  s 
CALL  CONJ 


CALL 
CALL 

CALL 
ARG2 

CALL 
CALL 
CALL 
CALL 
CALL 
DO  10 


CONJ 
CONJ 
XAR  ( 
s  ARG 
TAU  ( 
CONJ 
CONJ 
TAU  ( 
CONJ 
I  = 


CALL  TAU 

CALL  GUO 

CALL  XAB 

CALL  CONJ 


CONjGfALP(NN)  ) 
<4*NQ,U0(l.NPl)iV(l»NPl> > 
(4,U0»UO(l,NN>) 
<ALP(NN),UO(l»NN),TAO(l»NN)) 

(NN.N1.ALP.AA, IN,  I  I,KK) 

(NNPl, Nl, ALP) 

CONJG(ALP<NN) ) 
<N2iALP(NNPl)iALP(NNPl) ) 
(4*N2,U0<1» NNPl >#V(li NNPl) ) 
(4,V(l#NN)#U0(liNN) ) 

NNP1#1,N1#N#1) 

2(5) 

V(l',NN),TAO(l,N)  ) 

< 4*N2, V < li NNPl )iUO<1, NNPl)) 

(4, V(1,NN), V(1,N) ) 

V(l,N),TAO(l,NN) ) 

<4,U0(1,NN),V(1,NN) ) 

NNJP1.N1 

UO(l, I),TAO(l# I)) 

N2*1»U(1#NQ*1)#ALP(N) ) 

NP1,NN,NN,N1,3) 

<4,V(1,N),V(1.NN) ) 
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ARG2  =  ARC32(5> 

CALL  *AH  (NNP1  ,N,N1»M,2) 

RETURN 

END 


10 
20 


30 
40 
50 


60 
70 

80 
90 


SUBROUTINE  XAB  < N , M , NN , MM , KK ) 

THIS  ROUTINE  CAULS  THE  ROUTINES  WHICH  SOLVE  THE  DIFFERENCE 

EQUATIONS 

COMPLEX  UO,TAO,ARG,ARG2,U,T(2),ALP»BET,V 

COMMON  /A/  C(20  ), GAMMA, GST,ARG<5) ,ARG2<5) 

COMMON  /B/  UO(4. 1024) ,V<4, 1024), U(2,900>#TAO(2, 1024), D 

COMMON  /C/  ALP(1024),BET(1024) 

LL  =  1 

M2  =  M-l 

L  s  1*(KK-2)*(M*1-N) 

MMM  =  M  M ♦ 1 • M 

DO  70  I  s  N,NN 

IM1  s  1-1 

CALL  OARG  (V(l,  I  ).UO(l,M) ) 

ARG  =  ARG2 

GO  TO  (30,10,20),  KK 

IF  (I-M)  40,40,20 

CALL  VALU  (V(l, JMl),lO(l,M),U<l»L>» I , TAO ( 1, M ) , KK*2 ) 

LL  ■  2 

GO  TO  50 

CALL  GARG  ( UO < 1 , MM ) , V < 1 , I  Ml ) ) 

CALL  TAU  (UO(1,M),TAO(1,M) ) 

ARG2  s  ARG 

CALL  NEXT  (MMM,UO(l,M),TAO(l,M),U(l,D# I#M2.LL> 

IF  <KK,EQ.2.AND.LL.E0.2>  60,70 

BET( I )  s  U(2,MMM) 

CONTINUE 

GO  TO  (80,90,90),  KK 

CALL  GARG  (UO(l,MM), V(1,NN) ) 

RETURN 

END 


C 
C 


SUBROUTINE  VALU  ( U, V, X, I , U , K ) 

COMPUTES  X  AND  Y  FO»  THE  REGULAR  PART  OF  THE  SOLUTION  ALONG  THE 

SEE  CHARACTFRISTIC 

COMPLEX  U(4),V(4),X(2),W<4),T(2>,A,ALP,8ET,H,ES,B,B2 

COMMON  /C/  ALP(l024),aET(1024),B,B2,A,H 

COMMON  /D/  N(8) 


Qk 
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10 


20 


T  =  W(2) 

CALL  TAU  (V,W) 

T  s  ,5*<W(2)+T) 

H  s  ,5*ULP(  I  )-ALP(I-l)  ) 

A  s  ALP( I >»H 

IF  (K)  20,10,20 

CALL  JAK  (ALP(  I  >,X(2)»ALP) 

CALL  JAK  (A, BET. ALP) 

A  s  A*A 

X  x  X*ES(U(3),V(3),T)-2,*H*T*ALP 

RETURN 

T(2>  =  X(2) 

X(2)  =  BET( I ) 

A  =  A*A 

X  =  X+ES(U(3) ,V(3),T)-T*(X(2) 

RETURN 

END 


T(2)) 


10 


20 


30 
40 


SUBROUTINE  NEXT  ( N, U, Tl , XX , I  I , J J, LL > 

FINDS  THE  VALUES  NEEDED  FOR  THE  DIFFERENCE 

COMPLEX  U(4,2),X(2),T1 (2,1), 

DO  40  I  =  2.N 

IM1  =  1-1 

CALL  SOLVE  (Ud,  I),U(1,  IMD, 

CALL  TAU  (X,T2) 

DO  10  K  :  1,2 

T3(K)  =  ,5*{T1(K,  IMD+T2CK)) 

T2(K)  s  ,5*<Tl<KiI)*T2<K>> 

CALL  SOLVE  (U(l,  I  ),U(1, IMD 

CALL  SOLVE  ( U ( 3 ,  I  )  , U ( 3 , 

GO  TO  (30.20),  LL 

S  s  ES(U(3, I ),X,T2(2) ) 

CALL  SOLVE  (XX(1,  I  ),XX( 

CALL  SAVE  (U(l,  I  ),X,XX( 

U(4#  I  )  s  x(2) 

CALL  TAU  (U(l, 

U(3# I )  s  X 

RETURN 

END 


EQUATIONS 
T2(2),T3(2).ES,H,ARG,ARQ2»XX(2,1)»S,AL 


X»Tl(l,I),Tl(2,IMl),(t),,0.)) 


I  Ml  ) 


U( 
X» 


1,  D.T2 
-T2(2), 


T3(2) 

T3,  (0 


(0, 
.0. 


,0. 

)) 


) ) 


1# 
1, 


DiTIUi  I  )  ) 


1M1),XX(1, I >,-T2(2),-T3,S) 

I  ), I  I,  JJ+P 


SUBROUTINE  SAVE  (X,Y,Z, I, J) 

THIS  ROUTINE  PUTS  THE  SOLUTION  FOR  REAL  POINTS  OF  THE  t-OCOC-RAPH  ON 

TAPE  AND  SELECTIVELY  PRINTS  THE  SOLUTION 
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110 
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140 


COMPLEX  ALP»BET,Y(2),B,Z(2),B2,U0,U#TA0,V 

COMMON  /C/  ALP(1024),PET(1Q24),B.B2 

COM M0N/B/U0(4, 1024), V(4, 1024), U<2, 900), TAO(2»1024),D(600),A<4, 600) 

COMMON  /D/  N,M,  IN.MM.NRN,  I  I  , JJ,KK 

DATA  ER/,lE-5/ 

DIMENSION  X(2,2>,  AA(4) 

IF  (KK)  10,10,20 

IF  (I-J)  130,30,130 

IF  (  I  ,LT, I  I .OR. J,LT. JJ)  130,30 

M  =  M*l 

IF  (M-l)  80,40,80 

AA  s  X 

AA(2)  s  X(l,2) 

AA(4)  s  ( Y-R2)/ALP( I >*Z 

AA(3)  =  Y(2)/ALP( I )*Z<2> 

IF  (  I-I  I  )  50,60,50 

D  s  PSI  ( A,AA,D) 

DO  70  K  =  1,4 

A(K)  =  A  A ( K ) 

IF  (MOD( I+IN-I I. IN))  130,90,130 

A<i,M)  s  X 

A(2iM)  =  X(l,2) 

A(4#M)  c  (Y-B2)/ALP( I )*Z 

A(3#M)  s  Y(2)/ALP( I >*Z<2) 

D(M)  =  PSI  (  A(1,M-1)  ,  A(.1.,M)  ,D(M-1)  ) 

IF  (MODC  I  +  IN-I  I, IN) )  130,90,130 

IF  (MOD(M, IN)-1)  100,110,100 

IF  (J+IN-MM)  130,110,110 

PRINT  140,  (A(K,M),K  =  1,4),D(M) 

IF  (ABS(X(2,1) )+ABS(X(2#2) )-ER)  130,130,120 

PRINT  140,  X(2,1),X(2,2) 

RETURN 

FORMAT  (5F16.10) 

END 


10 
20 


SUBROUTINE  TAU  (X,Y) 

COMPUTES  LAMBDA*  AND  LAMBDA-  GIVEN 

COMMON  /A/  C(20), GAMMA, QST,ARG(5) 

COMPLEX  X(2),Y(2),Q,CS,ARG,AR 

O  =  X(1)*X(1)*X(2)*X<2) 

CS  =  QST-C(9)*Q 

Q  s  CSQRT(CS*(Q-CS) ) 

CALL  BRAN  (Q,ARG,AR,J) 

IF  (J)  10,20,10 

PRINT  30,  X,CS,Q, ARG(1),AR 

ARG  =  AR 


U  AND  V 
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30 


Y(2)  =  CS-X*X 


Y(l)  - 

Y(2>  = 

RETURN 

FORMAT 
END 


(X(1)*X(2)*Q)/Y<2) 
(X(1)#X(2)-Q)/Y(2) 

(4E2fl.l0) 


10 


20 

30 

40 
50 


SUBROUTINE  RRAN  <X,A,C.J) 

SELECTS  THE  CORRECT  BRANCH  OF  THE  SOLARE  ROOT 

COMPLEX  X,A,C,B 

DIMENSION  D<2) 

EQUIVALENCE  (  R  ,  I)  > 

DATA  E/,1/ 

J  =  0 

C  s  X/CA8S(X) 

B  s  C-A 

Y  s  D*D*D<2)*D<2>-1. 
IF  (Y)  20,30,10 

X  s  -X 
C  s  -C 

Y  =  -Y 

IF  <Y*E>  40,40,30 

J  =  2 

PRINT  50,  X,A,C 

RETURN 

FORMAT  (4E20.12) 

END 


DATA 


SUBROUTINE  JAK  (X.Z.ZP) 

COMPUTES  THE  INITIAL  CHARACTERISTIC 

COMMON  /A/  C<20>, GAMMA, OST 

COMMON  /E/  CC(4) 

COMPLEX  X,Y,Z,ZP,S 

S  s  C<ia)*X*Xt-C(l9) 

Y  x  CSQRT(S) 

Z  x  C(20)*(Y/X+C(18)*X/Y)-1,/X 
ZP  s  (1,-C(19)/(C(20)*S*Y))/(X*X) 
S  x  CC(4)-X 

Y  =  CC(?)»X 

Z  x  Z*CC*CL0G<Y/CC(2))*CC<3)*CL0G(S/CC(4)) 

ZP  «  ZP-CC/Y-CC(3)/S 

RETURN 

END 
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C 
C 


COMPLEX  FUNCTION  ES<L,V,T) 

COMPOTES  THF  I NHOMOGENEOUS  TERM  OF  THE  DIFFERENCE  EQUTICNS  FOR  X 

AMU  Y 

COMPLEX  U(2)fV(^),TiAl.,H,ALP#BET.B,R2 

COMMON  /C/  ALP( 1024) ,RET(1024)  ,R,R2, AL#H 

ES  =  H*(U*V-B*T*<U(2)+V(2)))/AL 

RETURN 
END 


SUBROUTINE  SOLVE  ( X. Y, Z  .  A.H.C > 

SOLVES  THE  SYSTEMS  OF  T^O  SIMULTANEOUS 

COMPLEX  S»T.A,B.C»X(2>,Y(2),Z(2> 

S  =  X(1)-A*X(2)*C 

T  =  Y(1)-R*Y<2) 

Z(2)  *  (T-S)/(A-B) 

Z<1>  =  S+A*Z(2) 

RETURN 

END 


LINEAR  ALGEBRAIC  EOLATIONS 


C 
C 


10 


SUBROUTINE  OUO  (N.U.A) 

STORES  THE  CHARACTERISTIC  INITIAL  DATA 

SOLUTION  ON  THE  ETA  CHARACTERISTIC 

COMPLEX  U(2,1),A,B,C 

CALL  JAK  (A.BfC) 

DO  10  I  =  1,N 

U(2i  I  )  -    B 

U<1*  I  )  s  (0..0,  ) 

RETURN 

END 


FOR  THE  REGUALR  PART  OF  THE 


10 


SUBROUTINE  QARG  (A,R) 

SETS  THE  FOUR  DIMENSIONAL 

COMPLEX  A(4)  tB(4) 

DO  10  I  =  1,4 

8(1)  =  A(  I  ) 

RETURN 


VECTOR  B  EGUAL  TO  A 
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END 


10 


SUBROUTINE  CONJ  (M,X,Y) 

SETS  THE  M-OIMENSIONAL  ARRAY  Y  EQUAL  TO  THE  COMPLEX  CONJLGATE  OF  X 

COMPLEX  X(1),Y(1  ) 

DO  10  1  =  1,M 

Y( I )  =  CONJG(X( I > ) 

RETURN 

END 


FUN 

COM 

DIM 

U  = 

V  = 

DY 

DX 

Q  s 

PSI 

RET 

END 


CTION  PSI  (X.Y.R) 

PUTES  THE  STREAM  FUNCTION  PSI 

ENSION  X(4),  Y<4> 

,5*(X+Y) 

l5*(X(2>*Y<2)  ) 
a  Y(4)-X<4) 
=  Y(3)-X(3) 

U*U*V*V 

s  R*RHO<Q)*(U*DY*V*DX) 
URN 


FUNCTION  RHO  (Q) 
COMPUTES  THE  DENSITY  RHO 
COMMON  /A/  C<20> 
RHO  =  (li=C(5)*Q)**C(10) 

RETURN 
END 


(EM, A,B) 

PAGE  CF  THE 
A(4) 


OUTPUT 


SUBROUTINE  TITLF 
PRINTS  THE  TITLE 
DIMENSION  EM(2), 
PRINT  10,  EM,B,A 
RETURN 
10  FORMAT  <lH,///2f)X,*  COMPRESSIBLE  FLOW 
1R*///25X*THF  MACH  NUMBER  AT  INFINITY  = 


ABOUT  AN  ELLIPTICAL  CYLINDE 
*#E10,3///30X*E  =  *jF5.2»8X 


2*  A  s  *,F5,2///27X4F8,2////10XlHU15X.lHV.15XtHXl5XlHY14X3l-PSI////> 
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END 


10 


WHICH  APPEAR  IN  THF  OLTPUT 


SUBROUTINE  HEAD 

LABELS  THE  COLUMNS 

PRINT  in 

RETURN 

FORMAT     UHl,9XlHU.i'!»XiHVl*XlHXlSXlHY,14X3HPSI////> 

END 


C 
C 
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20 

30 
40 
50 


PRO 
TH1 

COM 

COM 

DIM 

DRA 

REa 

REW 

REW 

REa 

PRI 

IF 

REA 

CAL 

CAL 

REA 

SFR 

XM 

SF 

CAL 

CAL 

CAL 

CAL 

CAL 

CAL 

CAL 

CAL 

PRI 

CAL 

CAL 

CAL 

FOR 

FOR 

FOR 


ORAM  SETGRF( INPUT, OU TPUT#TAPE1,TAPE2) 

S  PART  FINDS  THE  AIRFOIL  AND  DOES  THE  PLOTTING 

MON  /A/  0(22) 

MON  /D/  SF.SFR, YM, YY,NN 

ENSION  R(8) ,  TEXT(4) 

WS  AND  LABELS  AXES  AND  SETS  PLOT  PARAMETERS 

DS  THE  DATA  TO  BE  PLOTTED 

IND  1 

IND  2 

U  (1)  C,NNi IN 

NT  30,  C(17),C(18) 

(NN)  10,20,10 

D  6  0  »  R 

L  PLOTS  (2,1) 

L  PLOT  (0.  ,-11. ,-3) 

D  40,  YM,XM 

s  6,/YM 
=  SFR*(2.*XM-YM)/2, 
=  1,/SFR 

PAXIS  (2. ,5. ,5, ) 

XAXIS  (2. ,4,5,7, ) 

SYMBOL  (3. ,1.5, ,14,R,Q,,3Q) 

SYMBOL  (2, . .75,  ,14,R(4),0 ,,50) 

PLOT  (5.5+XM.4.5.-3) 

GRF  (IN.MN) 

PLOT(6.5-XM,0.,-3) 

PLOT  (0.,0.,99<?) 
50 


L 

L 

L 

L 

L 

L 

L 

L 

NT 

L  EXIT 

L  GRF  (  IM,NN) 

L  EXIT 

MAT  (30X.2F10 

MAT  (2F10.2) 

MAT(1H, ) 


4/) 


90 


p  a  rc b   16 


60 


FORMAT 
END 


(8A10) 


10 

20 
30 


40 

50 
60 
70 

80 
90 
100 
110 
120 
130 


140 


150 


160 

170 
180 

190 

200 


SUBROUTINE  GRF  MN.NN) 

FINDS  AND  PLOTS  POINTS  ON  THE  BODY 

PLOTS  THE  PRESSURE  DISTRIBUTION 

DIMENSION  A(4,600),  D(600),  CP(500), 

COMMON  /D/  SF,SS,YM 

DATA  EP,SR,K,LL/.1E-4,3,5,0,1/ 

READ  (1)  KK.Li ( (A(I,  J),  I  =  1,4), 

IF  (ENDFILE  1)200,20 

IF  (KK)  200,30,60 

K  s  -1 

XX  a  D(L) 

JF  <ABS(A(2,L>) 


X<500),  AA(5) 


D( J) , J  »  1,L) 


•  EP)  40,50,50 


ZO 
GO 
IF 
IF 
K 


XX-.1E-10 


TO  130 

(ABS<A(2) )-EP> 
(K)  70,80,100 
s  0 


120,130,130 


D( J)*XX 
.1E-12 


GO  TO  90 
READ  230,  XX 

XX  a  XX-D(L) 
DO  110  J  s  1, 

D(J) 

ZO  a 

DO  180  M  =  3,L 

M2  s  M-2 

CALL  INTERP  (  D < M2 ) , 70 » A ( 1 , M2 ) , AA , J, 4  ) 

GO  TO  (180,140),  J 

ANQ  =  180, *ATAN2(AA(2), AA)/3, 14159 

AA(9)  s  PE(AA) 

PRINT  240,  ( AA(N) ,N  s 

IF  (NN)  150,180,150 

X(Ll)  s  AA(3> 

CP(LD  b  -AA(5) 

CALL  SYMBOL  (SS*AA(3) 

IF  (KK)  160,170,160 

LL  *  LL*1 

GO  TO  180 

CALL  SYMBOL 

CONTINUE 

IF  (KK)  200 

CALL  BE  (K, 

GO  TO  10 

IF(LL,LE,2, 


1,5) 


SS*AA(4),  ,07,15,ANG,-1.) 


(SS*X(LL),SR+2,*CP(LL),.07,3,0,,-1) 

,10,390 
L,  A,  IN) 

O.NN.EQ.O)  220,210 
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210  CALL  PLHT  <0..SR,-3) 

X(|_L>  =  CP(LL)  =  0. 

CP(LL*1)  =  .5 

X(LL*1>  =  sr 

CALL  LINE  (y,CP,LL-l»l»C',3) 
220  RETURN 
230  FORMAT  (F12.8) 
240  FORMAT  (5F1-S.10) 

END 
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60 


SUBROUTINE  ME  ( J,N, A, IN> 

PLOTS  CHARACTERISTICS  IN  THE  SUPERSONIC 

DIMENSION  AA(2,400,27), A<4.1) 

COMMON  /D/  SF.SS 

J  =  J*l 

L  s  0 

DO  10  I  =  1,N, IN 

L  =  L*l 

AA(1, J,L>  =  A(3.  I) 

AA(2, J,L)  =  A(4,  I ) 

IF  (N-l)  50,50,20 

IF  (MOU( J, IN ) -1 )  40,30,40 

A(3»N*1)  =  A(4,.M*1)  =  0. 

A(3#N+2)  =  A(4,N+2)  s  SF 

CALL  LINE  ( A(3)  , A(4) ,N,4,0,3) 

RETURN 

M  =  J 

MM  =  IN*M 

L  =  M/IN 

DO  60  I  =  1,L 

K  =  MM- J*1N 

A  A  ( 1 ,  K  + 1  #  I  ) 

AA(1,K*2,  I  ) 


REGION 


AA(2,K  +  1,  I  ) 
AA(2»K*2, I > 


0. 
SF 


CALL  LINE 

RETURN 

END 


(  AAd.l,  I  ),AA(2,1,  I),K,2,0,3) 


FUNCTION  PE  (X) 

COMPUTES  THF  COEFFICIENT  OF  PRESSURE 

DIMENSION  X(2) 

COMMON  /A/  C(20  ),  GAMMA,  'JST 

Q  =  X*X+X(2)*X(?) 

PE  =  C*(  <QST-C(9)*Q>*RH'J(0)-C<4)  ) 
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RETURN 
END 


SUBROUTINE  INTERP  ( V, T, U, w, J, N ) 


10 
20 


30 
40 


COMPUTES  POINTS  ON  THE  BODY 
DIMENSION  Y(l),  U(l),  V<6), 


BY  QUADRATIC 
W(l) 


INTERPOLATION 


H(A,B,C,D,E) 

J  =  1 

XB 

XC 

IF 

IF 

J 

EA  = 
EB  = 
EC  s 
N2  = 
XA  s 
DO  30 
W(  I  ) 


(A*C-B*C)/E 


=  Y(2)«T 
■  Y(3>=>T 

CXB*XC)  20,10,40 
< (X*Z)  .EO.O.  )  40,20 
2 

=  Y(3)=Y(2) 
=  Y(3)=Y 
=  Y(2)=Y 

s  2*N 
Y-T 

I  =  1,N 
=  H(U(N*I  ),U(N2*I ) ,XC,XB,EA) 


V(  I>  - 
W(  I  )  = 

RETURN 
END 


H(U( I ),U(N2*J ),XC,XA,EB) 
H(V( I ),W(I),XB,XA,EC) 


FUNCTION  RHO  (O) 
COMPUTES  THE  DENSITY  RHO 
COMMON  /A/  C(20) 
RHO  s  (l,=C(5)*O)**C(10) 

RETURN 
END 


SUBROUTINE  XAXIS  (X,Y,S) 

DRAWS  THE  CHORD  AXIS 

J  =  S*l 

CALL  PLOT  (x,Y,3) 

XX  s  X-l 

CALL  PLOT  (X*S,Y,2> 
DO  10  I  =  1,J 
10  CALL  SYMBOL  ( XX*  I  . Y,  ,  14 , 16, 0 . , -1 ) 
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XX    =    XX  +  S+,<S 

CA|_L    SYMBOL 
RETURN 

END 


(XX. Y-,5,  ,14,1HX,0,,1) 


10 


SUBROUTINE  PAXIS  t X , Y , S  ) 

DRA^S  THE  COEFFICIENT  QF 

J  =  S+l, 

SP  =  S*1,*Y 

XX  =  X-.B 

<X,Y,3> 
(X, Y*S,P) 
1,  J 


PRESSURE  AXIS 


CALL  PLOT 
CALL  PLOT 
DO  10  I  = 
YY  =  SP-I 
VAL  =  ,5*< 1-3) 


CALL 
CALL 
CALL 
CALL 


SYMBOL 
NUMBER 
SYMBOL 
SYMBOL 


(X.YY, .14.15,0,1-1) 

(XX, YY. ,14,VAL,0..1) 
(X+.3.Y*S-,5, ,14,1HC»0 ., 
(X+.5.Y+S-.6, ,14,1HP,0,, 


1) 
1) 


RETURN 
END 
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COMPRESSIBLE  FLUW  ABCUT  AN  ELLIPTICAL  CYLINDER 


THE  maCH  NUMBER  AT  INFINITY  =   7.906E-01 


E  s    .50 


A  a   -  ,  4(1 


+  .40 


.59     +  .15    -1.50 
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PSI 


,7841014584 
.7798043333 

.78157444192 
.6024900546 

.6308799409 
.6735921987 
.852828i?!3o 
.6321303180 
.8121013821 


U. 0000000 n 00 
. 0469325*31 
.10U5/71989 
.l55bi20l84 

.2146/24321 
.2795805323 
.34085(16823 
,40*2/38635 
.4859078011 


-3.5314U65218 

-3.4320624301 

-3.2539960991 
-3.  049*007366 
-2.8595740236 
-2.  /064605252 
-2. 68493l4io3 
-2.O8U607509 
-2.6887196846 


0.0  000000000 
.3103392567 

.5103567593 
.5853598488 
.5532895157 

.4448328211 
.?9Q3Q99620 

.1774074603 

. J960939105 


0.0000000000 
.2009806441 

.3206190558 
,34826993i9 
.2982690278 
.1949175136 
.0853381455 
.0100435375 
-.0390999385 


OCT  2" 
Date  Due 


MAfc  2  &  1971 


2-jll^- 


ia|UCL^A. 


i: 


the 
of 
ort, 


ted 
itely 


to 
from 


behalf 
or 

of 
em- 
ir 
i- 
'or- 

itract 

with 


96 


c.l 


Korn 


Computation  of  shock-free 
transonic  flows  for  airfoil 
design. 


c.l 


Kern 


AUTHOR 

Computati 

T,-t£Fa,n sonic   f] 
foil   deL' 


MAR-e-srnfifj 


nj^^ti^rr 


N.Y.U.  Courant  Institute  of 
Mathematical  Sciences 

251  Mercer  St. 
New  York,  N.  Y.    10012 


